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Report  Title 

Mathematical  modeling  of  flow  through  vegetated  regions 

ABSTRACT 

Understanding  flow  processes  of  sea  and  fresh  water  through  complex 
coastal  regions  is  of  utmost  importance  for  a  number  of  applications  of  interest 
to  the  scientific  and  engineering  community,  including  wetland  health 
and  restoration,  inland  flooding  due  to  tropical  storms  and  hurricanes,  and 
navigation  through  coastal  waters.  In  such  regions,  the  existence  of  vegetation 
increases  flow  resistance,  which  is  a  major  factor  in  determining  velocity  and 
water  level  distribution  in  wetlands  and  inland.  Commonly,  the  momentum 
loss  due  to  vegetation  is  included  in  a  bottom  friction  term  in  the  model  equations; 
however,  such  models  may  oversimplify  the  complex  resistance  characteristics 
of  such  a  system.  With  recent  increases  in  computational  capabilities, 
it  is  now  feasible  to  develop  and  implement  more  intricate  resistance  models 
that  more  accurately  capture  these  characteristics. 

We  present  two  methods  for  modeling  flow  through  vegetated  regions. 

With  the  first  method,  we  employ  mathematical  and  computational  upscaling 
viii 

techniques  from  the  study  of  subsurface  flow  to  parametrize  drag  in  a  complex 

heterogeneous  region.  These  parameterizations  vary  greatly  depending 

on  Reynolds  number.  For  the  coastal  flows  in  which  we  are  interested  the 

Reynolds  number  at  different  locations  in  the  domain  may  vary  from  0(1) 

to  0(1000),  so  we  must  consider  laminar  and  fully  turbulent  flows.  Large 

eddy  simulation  (LES)  is  used  to  model  the  effects  of  turbulence.  The  geometry 

of  a  periodic  cell  of  vegetative  obstacles  is  completely  resolved  in  the 

fluid  mesh  with  a  standard  no-slip  boundary  condition  imposed  on  the  fiuidvegetation 

boundaries.  The  corresponding  drag  coefficient  is  calculated  and 

upscaling  laws  from  the  study  of  inertial  flow  through  porous  media  are  used 

to  parametrize  the  drag  coefficient  over  a  large  range  of  Reynolds  numbers. 

Simulations  are  performed  using  a  locally  conservative,  stabilized  continuous 

Galerkin  finite  element  method  on  highly-resolved,  unstructured  2D  and  3D 

meshes. 

The  second  method  we  present  is  an  immersed  structure  approach.  In 
this  method,  separate  meshes  are  used  for  the  fluid  domain  and  vegetative 
obstacles.  Taking  techniques  from  immersed  boundary  finite  element  methods, 
the  effects  of  the  fluid  on  the  vegetative  structures  and  vice  versa  are 
calculated  using  integral  transforms.  This  method  allows  us  to  model  flow 
over  much  larger  scales  and  containing  much  more  complicated  obstacle  geometry. 
Using  a  simple  elastic  structure  model  we  can  incorporate  bending 
and  moving  obstacles  which  would  be  extremely  computationally  expensive 
for  the  first  method.  We  model  flexible  vegetation  as  thin,  elastic,  inextenix 
sible  cantilever  beams.  We  present  two  numerical  methods  for  modeling  the 
beam  motion  and  analyze  their  computational  expense,  stability,  and  accuracy. 

Using  the  immersed  structure  approach,  a  fully  coupled  steady-state 
fluid-vegetation  interaction  model  is  developed  as  well  as  a  dynamic  interaction 
model  assuming  dynamic  fluid  flow  and  quasi-static  beam  bending.  This 
method  is  verified  using  channel  flow  and  wave  tank  test  problems.  We  calculate 
the  bulk  drag  coefficient  in  these  flow  scenarios  and  analyze  their  trends 
with  changing  model  parameters  including  stem  population  density  and  flow 
Reynolds  number.  These  results  are  compared  to  well-respected  experimental 
results.  We  model  real-life  beds  of  Spartina  altemiflora  grass  with  representative 
beds  of  flexible  beams  and  perform  similar  comparisons. 
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Understanding  flow  processes  of  sea  and  fresh  water  through  complex 
coastal  regions  is  of  utmost  importance  for  a  number  of  applications  of  in¬ 
terest  to  the  scientihc  and  engineering  community,  including  wetland  health 
and  restoration,  inland  flooding  due  to  tropical  storms  and  hurricanes,  and 
navigation  through  coastal  waters.  In  such  regions,  the  existence  of  vegetation 
increases  flow  resistance,  which  is  a  major  factor  in  determining  velocity  and 
water  level  distribution  in  wetlands  and  inland.  Commonly,  the  momentum 
loss  due  to  vegetation  is  included  in  a  bottom  friction  term  in  the  model  equa¬ 
tions;  however,  such  models  may  oversimplify  the  complex  resistance  charac¬ 
teristics  of  such  a  system.  With  recent  increases  in  computational  capabilities, 
it  is  now  feasible  to  develop  and  implement  more  intricate  resistance  models 
that  more  accurately  capture  these  characteristics. 

We  present  two  methods  for  modeling  flow  through  vegetated  regions. 
With  the  hrst  method,  we  employ  mathematical  and  computational  upscaling 
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techniques  from  the  study  of  subsurface  flow  to  parametrize  drag  in  a  com¬ 
plex  heterogeneous  region.  These  parameterizations  vary  greatly  depending 
on  Reynolds  number.  For  the  coastal  flows  in  which  we  are  interested  the 
Reynolds  number  at  different  locations  in  the  domain  may  vary  from  0(1) 
to  (9(1000),  so  we  must  consider  laminar  and  fully  turbulent  flows.  Large 
eddy  simulation  (LES)  is  used  to  model  the  effects  of  turbulence.  The  ge¬ 
ometry  of  a  periodic  cell  of  vegetative  obstacles  is  completely  resolved  in  the 
fluid  mesh  with  a  standard  no-slip  boundary  condition  imposed  on  the  fluid- 
vegetation  boundaries.  The  corresponding  drag  coefficient  is  calculated  and 
upscaling  laws  from  the  study  of  inertial  flow  through  porous  media  are  used 
to  parametrize  the  drag  coefficient  over  a  large  range  of  Reynolds  numbers. 
Simulations  are  performed  using  a  locally  conservative,  stabilized  continuous 
Galerkin  hnite  element  method  on  highly-resolved,  unstructured  2D  and  3D 
meshes. 

The  second  method  we  present  is  an  immersed  structure  approach.  In 
this  method,  separate  meshes  are  used  for  the  fluid  domain  and  vegetative 
obstacles.  Taking  techniques  from  immersed  boundary  finite  element  meth¬ 
ods,  the  effects  of  the  fluid  on  the  vegetative  structures  and  vice  versa  are 
calculated  using  integral  transforms.  This  method  allows  us  to  model  flow 
over  much  larger  scales  and  containing  much  more  complicated  obstacle  ge¬ 
ometry.  Using  a  simple  elastic  structure  model  we  can  incorporate  bending 
and  moving  obstacles  which  would  be  extremely  computationally  expensive 
for  the  first  method.  We  model  flexible  vegetation  as  thin,  elastic,  inexten- 
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sible  cantilever  beams.  We  present  two  numerical  methods  for  modeling  the 
beam  motion  and  analyze  their  computational  expense,  stability,  and  accu¬ 
racy.  Using  the  immersed  structure  approach,  a  fully  coupled  steady-state 
fluid-vegetation  interaction  model  is  developed  as  well  as  a  dynamic  interac¬ 
tion  model  assuming  dynamic  fluid  flow  and  quasi-static  beam  bending.  This 
method  is  verihed  using  channel  flow  and  wave  tank  test  problems.  We  calcu¬ 
late  the  bulk  drag  coefficient  in  these  flow  scenarios  and  analyze  their  trends 
with  changing  model  parameters  including  stem  population  density  and  flow 
Reynolds  number.  These  results  are  compared  to  well-respected  experimental 
results.  We  model  real-life  beds  of  Spartina  alterniflora  grass  with  represen¬ 
tative  beds  of  flexible  beams  and  perform  similar  comparisons. 
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Chapter  1 
Introduction 

1.1  Purpose  of  Study 

Understanding  flow  processes  of  sea  and  fresh  water  through  complex 
coastal  regions  is  important  for  a  number  of  applications  of  interest  to  the 
scientific  and  engineering  community;  e.g.,  wetland  protection  and  restoration 
projects,  navigation  through  coastal  waters,  and  projects  aimed  at  reducing 
the  effects  of  storm  surge  on  low-lying  coastal  areas.  Coastal  flow  processes  and 
navigation  are  complicated  by  the  presence  of  smaller-scale  features,  including 
stationary  and  moving  man-made  structures,  vegetated  marshlands,  barrier 
islands,  etc.,  and  the  interaction  of  these  features  with  the  larger  flow  domain 
which  may  include  the  continental  shelf  and  ocean  basins. 

In  the  applications  mentioned  above,  particularly  coastal  inundation 
due  to  storm  surge,  the  ability  of  these  smaller-scale  features  to  provide  re¬ 
sistance  to  flow  is  of  utmost  concern.  Current  mathematical  models  and 
parametrizations  of  flow  resistance  characteristics  due  to  flow  over  and  around 
porous  structures  or  through  heavily  vegetated  regions  are  inadequate  to  ac¬ 
curately  capture  these  physical  processes,  and  as  a  result,  numerical  models 
based  on  these  mathematical  representations  may  be  inaccurate.  Furthermore, 
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while  the  use  of  highly  resolved,  three-dimensional  Navier-Stokes  models  may 
avoid  some  of  the  parametrization  issues,  the  implementation  of  such  models 
over  large  coastal  regions  is  prohibitive,  and  in  many  cases,  unnecessary  to 
adequately  describe  the  flow  physics.  An  alternative  is  to  appropriately  cou¬ 
ple  mathematical/numerical  models  that  accurately  and  efficiently  capture  the 
primary  flow  characteristics  of  a  given  region  within  the  larger  flow  domain. 

In  this  thesis  we  focus  on  the  flow  of  water  through  a  vegetated  region. 
The  existence  of  vegetation  affects  the  flow  resistance,  which  is  a  major  factor 
in  determining  velocity  and  pressure  distribution.  The  main  impact  of  vege¬ 
tation  on  flow  is  form  drag,  which  results  in  momentum  losses  as  discussed  by 
Batist  et  al.  [4] .  Since  modeling  the  effects  on  flow  from  each  individual  plant 
is  impractical  in  a  large-scale  simulation,  where  computational  resolution  is  on 
the  order  of  tens  of  meters,  we  focus  on  methods  for  describing  and  capturing 
this  sub-grid  scale  phenomena  through  the  computation  of  drag  coefficients 
which  may  be  used  in  a  macro-scale  model. 

1.2  Background 

A  great  deal  of  experimental  and  computational  effort  has  been  focused 
on  characterizing  drag  effects  of  flow  through  vegetated  regions.  Much  exper¬ 
imental  work  has  focused  on  calculating  Manning’s  n  for  flow  over  different 
types  of  plants.  The  Manning  formula  is  a  highly-simplified  empirical  model 
for  fully  turbulent  open  channel  and  free-surface  flow.  Manning’s  n  incor¬ 
porates  bottom  friction,  bottom  shape,  vegetative  drag,  and  other  resistance 
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characteristics  into  one  value.  According  to  experimental  work  done  by  Ree 
and  Palmer  [87]  and  Petryk  and  Bosmajian  [79],  it  varies  greatly  depending  on 
properties  of  vegetation  and  is  highly  dependent  on  water  depth;  although,  this 
factor  is  commonly  ignored.  Thus  many,  including  Kandlec  [43]  and  Turner 
and  Chanmeesri  [105]  have  demonstrated  that  Manning’s  equation  is  inade¬ 
quate  for  describing  shallow  water  flow  through  vegetation  over  a  large  range 
of  velocities  and  for  a  variety  of  water  depths.  They  suggest  that  the  appro¬ 
priate  drag  models  would  best  be  site-specihc  functions  of  water  depth  and 
plant  shape,  distribution,  and  flexibility. 

In  this  report,  we  use  methods  from  the  study  of  subsurface  flow  in 
porous  media  and  high-resolution  modeling  of  flow  through  vegetated  regions 
to  calculate  upscaled  drag  coefficients,  which  vary  greatly  with  Reynolds  num¬ 
ber  and  water  depth,  that  may  be  used  in  large-scale  resistance  modeling  of 
flow  through  vegetated  areas.  Upscaling  is  the  process  of  analyzing  micro-scale 
dynamics  to  describe  larger  scale  phenomena.  Vegetated  domains  have  sim¬ 
ilar  geometry  to  subsurface  porous  media  domains.  We  use  simple  domains 
of  arrays  of  cones  and  cylinders  to  model  vegetation  structures.  Especially  at 
lower  Reynolds  numbers,  much  experimental  and  theoretical  research  has  been 
performed  for  upscaling  subsurface  flows.  Specihcally,  there  are  several  math¬ 
ematical  homogenization  methods  for  upscaling  viscous  flow  in  porous  media. 
Homogenization  is  the  method  of  mathematically  deriving  effective  equations 
to  describe  dynamics  at  larger  scales  by  satisfying  to  some  degree  the  dynam¬ 
ics  at  the  smaller  scales.  For  higher  Reynolds  numbers,  inertial  effects  and 
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turbulence  are  important  factors  in  flow  characteristics.  Reynolds  Averaged 
Navier-Stokes  (RANS)  models  and  modern  high-resolution  techniques  such  as 
large  eddy  simulation  (LES)  and  direct  numerical  simulation  (DNS)  are  useful 
in  resolving  these  factors. 

The  most  basic  model  of  flow  through  porous  media  is  Darcy’s  law, 
determined  experimentally  by  Henry  Darcy  in  1856.  Darcy’s  law  is  only  valid 
for  slow,  viscous  flow,  which  is  sufficient  for  most  subsurface  modeling  scenar¬ 
ios.  Darcy’s  law  contains  a  hydraulic  conductivity  tensor,  which  measures  the 
ability  of  the  media  to  permit  flow,  which  must  be  known  before  the  law  can  be 
used  for  predictive  simulations.  There  are  both  experimental  and  theoretical 
approaches  to  calculating  the  hydraulic  conductivity  of  a  porous  medium.  One 
method  is  to  measure  velocity  and  hydraulic  gradient  directly  and  use  Darcy’s 
law  to  calculate  the  hydraulic  conductivity.  It  is  also  common  for  hydraulic 
conductivity  to  be  estimated  using  empirically  derived  formulas.  Descriptions 
and  analyses  of  these  methods  are  done  by  Bear  [5].  Darcy’s  law  is  theoreti¬ 
cally  derived  by  mathematical  homogenization  of  Stokes  flow.  The  process  of 
homogenization  explained  by  Hornung  [37]  describes  a  theoretical  method  of 
determining  hydraulic  conductivities.  A  more  recent  in-depth  mathematical 
treatment  of  homogenization  theory  has  been  done  by  Mikelic  [70]. 

For  flow  through  vegetated  domains  we  generally  cannot  assume  slow, 
viscous  flow,  so  Darcy’s  law  is  an  inadequate  model.  For  higher  velocity  flows 
where  inertial  forces  matter,  there  exist  “non-Darcy”  models  which  include 
higher  order  terms.  These  laws  suggest  a  nonlinear  relationship  between  fluid 
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velocity  and  hydraulic  gradient  at  higher  Reynolds  numbers.  Forchheimer 
[26]  was  likely  the  first  to  propose  a  non-Darcy  relationship.  The  Darcy- 
Forchheimer  Equation  is  Darcy’s  law  with  an  additional  quadratic  velocity 
term.  The  Darcy- Forchheimer  Equation  has  been  observed  experimentally  for 
moderate  Reynolds  numbers  by  many  researchers.  There  have  been  many  at¬ 
tempts  at  deriving  the  non-Darcy  laws  using  mathematical  homogenization; 
however,  these  studies  such  as  those  by  Lions  [58],  Sanchez-Palencia  [91],  Mei 
and  Auriault  [67],  Balhoff  et  al.  [2],  and  Rasoloarijaona  and  Auriault  [83]  show 
the  quadratic  terms  canceling,  giving  a  cubic  law.  This  cubic  law  has  been 
verihed  numerically  for  very  small  Reynolds  numbers  by  solving  the  Navier- 
Stokes  equations  by  Colaud  et  al.[14],  Rojas  and  Koplik  [89],  Balhoff  et  al. 
[2],  and  Koch  and  Ladd  [50].  However,  this  has  only  been  verihed  where 
Darcy’s  law  is  approximately  valid,  and  the  quadratic  law  seems  to  be  more 
appropriate  for  higher  Reynolds  number  hows.  Some  methods  provide  a  result 
with  quadratic  and  cubic  terms.  Some,  including  Scheidegger  [93]  have  also 
used  power  laws  to  ht  empirical  data.  Garibotti  and  Peszyhska  [30]  provides 
computational  analysis  of  non-Darcy  how  through  porous  media  using  a  h- 
nite  diherence  method.  There  are  also  “computational  upscaling”  techniques. 
These  homogenization  methods  are  not  based  on  mathematical  derivation,  but 
are  performed  by  analyzing  results  from  solving  the  full  Navier-Stokes  equa¬ 
tions  on  hne  grids  as  seen  in  work  by  Peszyhska  et  al.  [78],  Balhoh  et  al.  [2], 
and  Balhoh  and  Wheeler  [3]. 

Higher  Reynolds  number  how  characteristics  through  vegetated  do- 
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mains  has  been  studied  a  great  deal  with  regard  to  atmospheric  flow  through 
plant  canopies  and  flow  through  vegetated  open  channels.  Dawson  and  Charl¬ 
ton  [19]  give  an  extensive  bibliography  of  earlier  work  in  this  held,  and  Nepf 
et  ah  [73]  [57]  [75]  provide  excellent  in-depth  experimental  analysis  of  drag, 
turbulence,  and  diffusion  for  flow  through  arrays  of  cylinders  and  dense  beds 
of  grasses  at  high  Reynolds  numbers.  Much  of  the  numerical  modeling  effort 
has  focused  on  developing  turbulence  closure  schemes  to  the  governing  equa¬ 
tions.  There  are  two  main  approaches  to  providing  turbulence  closure.  One 
method  uses  the  Reynolds  Averaged  Navier-Stokes  (RANS)  equations  which 
requires  a  model  for  Reynolds  stresses  to  provide  turbulence  closure.  The 
other  method  uses  Large  Eddy  Simulation  (LES)  turbulence  models  to  solve 
the  Altered  Navier-Stokes  equations. 

RANS  methods  have  been  used  with  varying  degrees  of  success.  Chris¬ 
tensen  [13]  found  that  simple  mixing  length  closure  methods  can  be  successful 
models  for  flow  through  simple  domains.  However,  Wilson  and  Shaw  [111] 
acknowledge  that  first  order  RANS  closure  schemes,  while  simple,  do  not  pro¬ 
vide  results  that  adequately  match  empirical  data  through  more  complex  do¬ 
mains.  They  propose  a  higher  order  closure  scheme  for  a  spatially  as  well  as 
temporally  averaged  version  of  the  governing  equations,  resulting  in  a  one¬ 
dimensional  representation  of  the  problem.  Raupach  and  Shaw  [86]  extend 
this  work,  proposing  a  method  of  obtaining  momentum  and  energy  equations 
in  multi-connected  flows.  These  equations  capture  different  momentum  and 
dissipation  terms  resulting  from  the  three-dimensional  nature  of  this  flow. 
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Raupach  et  al.  [85]  validate  this  model  experimentally,  using  aluminum  strips 
to  model  vegetation. 

The  most  commonly  used  and  most  useful  RANS  closure  schemes  are 
two  equation  methods.  These  schemes  involve  solving  two  transport  equations 
for  turbulent  kinetic  energy  and  dissipation  whose  solutions  dehne  an  eddy 
viscosity.  These  equations  contain  empirically  calculated  constants.  Burke 
and  Stolzenbach  [11]  introduce  drag- related  source  terms  to  model  the  effect 
of  vegetation.  The  most  commonly  used  two  equation  closure  models  are  the 
k  —  e  and  k  —  ui  models.  Lopez  and  Garcia  [60],  [61]  give  a  full  analysis  of  these 
schemes  and  show  their  ability  to  predict  3D  flow  patterns.  Dehna  and  Bixio 
[21]  shows  that  a  k  —  e  model  may  accurately  predict  flow  patterns  and  eddy 
viscosities,  but  may  poorly  predict  more  complex  turbulence  characteristics. 

As  an  alternative  to  RANS  models,  LES  can  be  used  to  solve  the  hltered 
Navier-Stokes  equations  and  can  provide  an  almost  complete  description  of  the 
flow,  while  not  requiring  empirically  derived  transport  equations  to  be  solved. 
Early  simulations  of  flow  and  turbulent  structures  above  forests  by  Moeng 
[71],  Shaw  and  Schumann  [95],  Kanda  and  Hino  [44],  and  Dwyer  et  al.  [24] 
show  that  important  turbulence  structures  cannot  be  captured  using  a  RANS 
model,  but  can  be  captured  by  an  LES  model.  LES  of  channel  flow  through 
vegetated  channels  is  presented  by  Cui  and  Neary  [15],  Stoesser  et  al.  [98], [99] 
and  Palau  and  Stoesser  [76].  In-depth  analyses  of  turbulence  statistics  and 
temporally-averaged  characteristics  of  these  flows  demonstrate  the  superiority 
of  LES  to  RANS  in  capturing  hne-scale  flow  qualities.  However,  LES  may 
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require  significantly  higher  resolution  than  RANS. 

Plant  flexibility  can  be  an  important  factor  in  resistance  characteris¬ 
tics.  Kouwen  and  Unny  [51]  performed  experiments  that  simulate  vegetation 
with  plastic  strips.  They  propose  that  there  exist  three  basic  types  of  flow 
through  vegetation:  when  the  vegetation  remains  erect  and  stationary,  when 
the  vegetation  undergoes  a  waving  motion,  and  when  the  vegetation  is  com¬ 
pletely  bent  over.  Similar  results  were  found  by  Gourlay  [32]  using  real  grass 
instead  of  plastic  strips.  Huai  et  al.  [38]  show  that  waving  vegetation  influ¬ 
ences  stream-wise  velocity  and  Reynolds  stresses.  There  have  been  numerical 
studies  incorporating  plant  flexibility.  The  approach  in  these  studies,  proposed 
by  Kutija  and  Hong  [53],  is  to  assume  a  flexible  ID  stem  and  a  velocity  prohle 
and  apply  classical  cantilever  beam  theory  to  simulate  bending.  In  the  region 
below  the  “effective  vegetation  height”  a  standard  drag  force  is  applied  and  a 
RANS  or  LES  technique  is  used  to  model  the  flow.  2D  LES  simulation  of  this 
type  of  model  was  performed  by  Ikeda  et  al.  [40].  Erduran  and  Kutija  [25]  use 
a  3D  RANS  technique  with  a  combination  of  mixing  length  and  eddy  viscosity 
turbulence  closure  schemes.  They  also  propose  a  quasi-3D  coupling  with  the 
shallow  water  equations.  Li  and  Xie  [54],  building  on  work  involving  stiff  vege¬ 
tation  by  Li  and  Zeng  [55],  add  to  the  flexible  vegetation  model  by  assuming  a 
thin  plate  of  “foliage”  and  use  a  3D  LES  scheme  with  Smagorinksy  turbulence 
closure.  While  these  models  have  been  successful  for  meso-scale  modeling  of 
fully  turbulent  free  surface  flow  over  submerged  and  non-submerged  vegeta¬ 
tion,  they  do  not  consider  a  valid  drag  model  for  laminar  or  transitional  flow 
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or  micro-scale  flow  and  do  not  capture  the  drag  effects  due  to  the  complex 
geometries  of  vegetation  at  the  micro-scale. 

In  this  work,  we  approach  vegetative  resistance  in  terms  of  drag  and  a 
variable  drag  coefficient,  and  numerically  simulate  flow  through  vegetation  for 
a  wide  range  of  Reynolds  numbers  with  two  main  approaches.  The  hrst  ap¬ 
proach  involves  fully  resolving  the  geometry  of  the  porous  domain  containing 
vegetative  obstacles.  For  low  Reynolds  number  flows,  drag  can  be  modeled 
with  techniques  from  the  study  of  flow  through  porous  media  such  as  Darcy’s 
law  and  non-Darcy  models.  We  look  at  the  theory  behind  Darcy’s  law  and  sim¬ 
ple  mathematical  homogenization  theory,  including  a  computational  method 
of  calculating  hydraulic  conductivities.  We  present  the  Darcy-Forchheimer 
equation  and  other  non-Darcy  laws  and  discuss  their  derivations  and  experi¬ 
mental  bases.  In  higher  Reynolds  number  flows,  eddies  form  and  turbulence 
becomes  important.  The  drag  coefficients  behave  much  differently  in  this  flow 
than  in  the  low  Reynolds  number  region.  We  incorporate  LES  using  the  sta¬ 
bilized  dynamic  Smagorinsky  eddy  viscosity  model,  allowing  us  to  realistically 
consider  the  effects  of  turbulence.  Using  a  stabilized  hnite  element  method 
that  we  explain  in  depth  in  the  next  chapter,  we  solve  the  equations  of  motion 
on  complex,  vegetated  2D  and  3D  domains  with  highly-resolved  meshes  with 
periodic  boundary  conditions.  This  allows  us  to  perform  a  qualitative  analysis 
of  drag  over  a  large  range  of  Reynolds  numbers.  We  also  analyze  the  accuracy 
of  Darcy’s  law  and  non-Darcy  equations,  including  the  ranges  of  Reynolds 
numbers  over  which  they  are  valid.  An  analysis  of  drag  coefficients  at  high 
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Reynolds  numbers  is  also  performed.  We  analyze  the  drag  coefficient  and  drag 
effects  of  vegetative  resistance  over  a  range  of  Reynolds  numbers  from  slow, 
viscous  flow  to  fully  turbulent  flow. 

Our  second  approach  is  an  immersed  structure  method.  With  this 
method,  we  use  separate  meshes  for  the  fluid  and  vegetative  obstacles  and  use 
ideas  from  fluid-structure  interaction  modeling  to  map  the  effects  of  the  two 
on  each  other.  With  this  method,  we  can  model  rigid  and  flexible  vegetation. 
Plant  flexibility  is  another  important  factor  in  vegetative  drag.  We  model 
flexible  vegetation  as  elastic  inextensible  cantilever  beams.  We  present  the 
theory  behind  such  beam  models  and  some  numerical  techniques  for  solving 
them.  We  develop  a  method  for  coupling  the  fluid-vegetation  systems  and 
test  its  accuracy.  We  use  this  method  to  calculate  bulk  drag  coefficients  for 
channel  flows  and  wave  tanks  containing  vegetative  obstacles.  Replicating  the 
geometry  of  experimental  setups  as  closely  as  possible,  we  verify  the  presented 
methods  by  comparing  to  experimental  results. 


10 


31 


Chapter  2 
Flow  Model 


In  this  chapter,  we  develop  the  model  that  is  used  in  almost  all  of 
the  flow  simulations  in  this  dissertation.  Beginning  with  the  fundamental 
Navier-Stokes  equations,  we  explain  the  formulation  of  a  continuous  Galerkin 
finite  element  method  using  variational  multiscale  stabilization.  We  explain 
the  importance  of  using  a  proper  turbulence  model  and  present  the  dynamic 
Smagorinksy  that  we  prefer  to  use.  Finally,  we  verify  our  flow  model  with  a 
difficult  benchmark  problem. 

2.1  Introduction 

The  standard  model  for  incompressible  Newtonian  flow  is  the  Navier- 
Stokes  equations.  This  is  a  system  of  conservation  laws  for  conservation  of 
mass  and  momentum.  In  conservative  form  they  are 

V-v  =  0  (2.1) 

^  + V-  (v®v-z/(Vv  +  Vv^)) +  — -g  =  0  (2.2) 

OT  P 

where  v  is  the  velocity,  p  is  the  pressure,  p  is  the  fluid  density,  u  is  the  kinematic 
viscosity,  and  g  is  a  body  force  including  gravity. 
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2.2  Variational  Formulation 

Let  f2  be  the  domain  of  the  fluid  region.  Suppose  that  the  boundary, 
dfl  has  two  partitionings  {r^,r^}  and  Suppose  that  we  have  the 


boundary  conditions 

p  =  Pd  on  (2.3) 

V  ■  n  =  on  (2.4) 

V  =  Vo  on  (2.5) 

[v  (g)  V  —  z/(x)(Vv  +  Vv^)]  ■  n  =  h)(r  on  r)(r.  (2.6) 

Assume  that  we  have  initial  conditions 

p(x,0)  =  po  (2.7) 

v(x,0)  =  vq.  (2.8) 

Dehne  the  trial  function  spaces 

VP  =  {w\w  e  H\n)  andw  =  pd  onTPj}  (2.9) 

=  {w|w  G  and  w  =  vo  on  L^}.  (2-10) 


Now,  where  t  G  [0,T]  is  time,  we  define  the  time-dependent  trial  function 
spaces:  Vf  =  {w\w{t)  =  V^}  and  V)’  =  {w|w(t)  G  V^}.  Dehne  the  test 
function  spaces  =  {w\w  G  and  w  =  0  on  L^}  and  =  {w|-w  G 

and  w  =  vo  on  0}.  We  multiply  the  equations  in  (2.1)  by  the  respec¬ 
tive  test  functions. 

[  {W  ■  v)wPdn  =  0  (2.11) 
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[  f^  +  V-(v®v  -  z/(Vv  +  Vv^))  +  ^  -  g 
JQ  \<^  P 


w^dn  =  0.  (2.12) 


Integrating  by  parts,  we  have  the  variational  formulation:  hnd  p  &  Vf  and 


V  G  Vj  such  that 


-  /  V  ■  WvjPdn 

Jq 


hPdT  WwP  e  WP, 


pp 
^  N 


(2.13) 


/  -  V  0  V  ■  /  z/(x)(Vv  + Vv"^)  • 

Jn  ^  Jn 


g - Vp  w^d^d 

'n  \  P 


h^j^dT  yw^  e  W" 


(2.14) 


2.3  Continuous  Galerkin  Formulation 

Let  be  a  simplicial  mesh  in  R”  with  n  =  2,3  containing  N^.  ele¬ 
ments,  {fie}  foi'  e  =  1,2,  Let  hg  be  the  diameter  of  fie-  Let  -P^(fle)  be 

the  set  of  polynomials  of  order  k  on  fie-  Dehne  the  discrete  spaces: 


r,f  =  {Pk  e  r”  n c»(sj)|pi  €  Pi(n,)}  (2.15) 

K  =  {<  6  W'" n c"(S)|u>!;  €  P^.(SJ,)}  (2.16) 

Yl  =  {vt  6  V”  n  C»(S)|v»  €  P!).(n.)}  (2.17) 

w;;  =  {v»  e  w”  n  c"(n)|vft  e  p^.(si.)}.  (2.18) 

The  classical  continuous  Galerkin  formulation  is:  hnd  ph  G  V)f  and  G  V]( 
such  that 

-  /  Vft  ■  V<dG  =  -  f  hPdT,  G  WP  (2.19) 

Jq 
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dvh 

dt 


wr  -  v/j  0  v/j  ■  VwldQ  = 


/  z/(x)(Vv,j  +  Vv^)  ■ 

'n 

^  (^g  - 


h^dr,  Vw;;  G  wi  (2.20) 


2.4  Variational  Multiscale  Formulation 

When  the  diffusive  effects  are  much  smaller  than  advective  effects,  the 
classical  continuous  Galerkin  method  tends  to  be  oscillatory.  It  becomes  es¬ 
sential  to  add  stabilization  to  the  formulation.  The  dominant  technique  in 
computational  fluid  dynamics  is  to  perform  a  multiscale  decomposition  of  the 
solution  to  derive  a  stabilized  method.  Stabilized  methods  tend  to  reduce  the 
oscillations  in  the  solution  without  dramatically  adding  computational  cost. 
Common  multiscale  decomposition  approaches  include  the  multiscale  hnite 
element  method,  the  subgrid  upscaling  technique,  and  the  mortar  method. 
These  methods  incorporate  small-scale  heterogeneity  into  the  formulation. 

The  variational  multiscale  (VMS)  method,  developed  by  Hughes  [39], 
is  a  method  that  decomposes  the  problem  into  a  grid-scale  problem  and  a 
subgrid-scale  problem.  We  derive  the  VMS  method  for  a  system  of  nonlinear 
conservation  laws 

^  +  V  •  (f(u)  -  D(u)Vu)  =  0  (2.21) 

where  u(x)  is  the  solution  vector,  x  G  C,  and  t  G  (0,T].  We  call  f  the 
advective  flux  and  D  the  diffusive  flux.  Let  the  boundary  of  the  domain  dVt 


14 


35 


be  decomposed  into  T ^  and  T ^  where  the  essential  and  natural  boundary 
conditions  respectively  are  imposed,  i.e. 

u  =  u  on  and  (2.22) 

(f  -  DVu)  ■  n  =  F.  (2.23) 

We  assume  the  initial  condition 

u(x,  t  =  0)  =  Uo(x),  X  G  n.  (2.24) 

This  VMS  formulation  for  nonlinear  conservation  laws  using  an  algebraic  sub¬ 
grid  scale  model  (ASGS)  is  shown  by  Juanes  and  Patzek  in  [42]. 

2.4.1  Weak  Form 

Dehne  the  variational  spaces 

V  ;=  {v  G  W|v  =  u  on  r £)}and 
Vo  ;=  {v  G  W|v  =  0  on  To} 

where  W  is  the  appropriate  Sobolev  space  depending  on  D.  Dehne 

((9tu,  v)  =  /  ^  ■  vdD,  (2.25) 

a(u,  v;w)  =  —  /  f(-w)  ■  VvdD -|-  /  D(w)Vu  ■  VvdD,  (2.26) 

Jn  Jn 

l(v)  =  -  [  F-ndT.  (2.27) 

The  weak  form  of  the  problem  for  each  t  G  (0,  T]  is  to  hnd  u  G  V  such  that 

((9tu,  v) -ha(u,  v;u)  = /(v),  Vv  G  Vq.  (2.28) 
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2.4.2  Galerkin  Method 

Let  \h  C  V  and  C  Vq  be  conforming  finite-dimensional  spaces  of 
piecewise  polynomials  on  a  finite  element  mesh.  Galerkin’s  method  is  to  hnd 
snch  that  for  each  t 

(dju/,,  V,,) -k  a(ufe,  Vfe;u,i)  = /(vfe),  Vv/^  e  V,i,o-  (2.29) 

2.4.3  Multiscale  Split 

Snppose  that  V  can  be  expressed  as  the  direct  snm  of  two  spaces 

V  =  V;,©  V.  (2.30) 

We  call  y h  the  space  of  resolved  scales  and  V  the  space  of  subgrid-scales.  Note 
that  V  is  an  inhnite  dimensional  space.  There  is  a  unique  multiscale  split 

u  =  u/i  +  u.  (2.31) 

For  nonlinear  problems,  Picard  or  Newton  iterative  methods  are  generally 
used.  For  these  methods  for  each  iterative  step  {k)  it  is  natural  to  represent 
the  solution  as 

uW  =  (2.32) 

Using  the  multiscale  split,  we  have 

(2.33) 

=  (2.34) 
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Assuming  we  have 

^  +  (2.35) 

For  simplicity,  we  drop  the  superscripts  u  ^  +  (5u. 

We  split  the  solution  and  the  test  function  into  the  multiscale  elements: 

u  =  u/j  +  (5u  G  V  =  V/i  ©  V  (2.36) 

V  =  v/j  +  V  G  Vq  =  V/i^o  ©  (2.37) 

(2.29)  can  now  be  split  into  a  grid  scale  problem 

{dt{uh  + 5\i),Vh)  =  a{uh  + 5\i,Vh;u  + 5u)  =  l{vh),  Vv/^  G  V,i,o  (2.38) 

and  a  subgrid  scale  problem 

((9t(u,i  +  (5u),v)  =a(u,i  +  (5u,v;u  +  (5u)  =  /(v),  Vv  G  V/^^o-  (2.39) 

2.4.4  Subgrid  Scale  Problem 

Consider  the  flux  term  in  (2.39)  as  a  sum  of  integrals  over  each  element. 
Integrate  by  parts  on  each  element  : 

a(u;i  +  5u,  v;  u  +  (5u) 

=  —  f  (f(u  +  (5u)  —  D(u  +  (5u)V  (u/i  +  (5u)))  ■  Vvdff 

g  J  lie 

=  —  f  V  •  (f(u  +  (5u)  —  D(u  +  (5u)V  {ufi  +  (5u)))  ■  vdhl 

e  Joe 

-  ^  [  ((f(u  +  5u)  -  D(u  +  (5u)V  (u/i  +  ((5u)))  ■  n)  ■  vdF 
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Assuming  continuity  of  flux  over  interior  element  boundaries 

—  f  ((f(u  +  (5u)  —  D(u  +  (5u)V  (uh  +  (5u)))  ■  n)  ■  vcfF 


^  /  F  ■  vdr  =  /(v) 

Approximate  the  total  flux  with  a  Taylor  expansion  about 


(2.40) 


f(u  +  (5u)  -  D(u  +  5u)V  (ufe  +  ((5u))  =  f(ufe)  -  D(u,i)Vu,i 


f'(u,i)5u  -  (D'(u,i)(5u)Vufe 

D(u;,)V(5fi)  +  (^(IH')-(2.41) 


We  define  the  linearized  advection-diffusion  operator 


:=  V  ■  [f  (u;,)v  -  (D'(u;,)v)Vu;,  -  D(u;,)Vv]  .  (2.42) 


We  can  also  write  this  as 


■■=  V  ■  [A(ufe)v  -  D(u/i)Vv] 


(2.43) 


where 


dfiiuh)  s^dDikiuh)^ 


(2.44) 


Now,  the  flux  term  in  the  subgrid  scale  problem  can  be  written  as 

a(u/j  +  (5u,  v;  u  +  (5u) 

~  f  V  •  (f(u  +  5u)  —  D(u  +  (5u)V  (u/j  +  ((5u)))  ■  vdhl 

e  Joe 

+  f  ■  ird^l  —  [  F  ■  vdr 


(2.45) 
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Assume  quasi-static  subscales: 


dtSu  ^  0.  (2.46) 

Define  the  grid-scale  residual 

^{uh)  :=  -dtUh  -  V  ■  (f(u;,)  -  D(uft)Vuft.  (2.47) 

The  subgrid-scale  problem  is  now 

f  JtfuhSu  ■  vdD  =  f  M{uh)  ■  vdD  Vv  G  V.  (2.48) 

We  cannot  expect  to  find  an  analytic  solution  to  this  infinite  dimensional 
problem.  We  use  an  algebraic  subgrid  scale  model  (ASGS)  for  an  algebraic 
approximation  for  the  subscales: 

5u  ^  Tu^^iuh).  (2.49) 

For  a  one-dimensional  system,  using  linear  elements  and  an  element  diameter 
of  h 

=  (^^D(u;,)  +  ^A(u;,)^  (2.50) 

is  an  appropriate  ASGS  model. 

2.4.5  Grid  Scale  Model 

Define  the  adjoint  of 

:=  -A'^(u,)Vv  -  V  •  (D'^(u,)Vv)  (2.51) 
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and  the  associated  boundary  operator 

b:,v  :=  (D’'(uk)Vv)  .  n,  (2.52) 

Now  (2.38)  becomes 

{^tUh,Vh)  +  a{uh,^rh;uh)  +  '^  SudQ  +  '^  hl^Vh- SudT  =  l(vh) 

g  ^  g  ^Te 

(2.53) 

Vv/i  G  'Vhfi-  Notice  that  this  is  (2.29)  with  two  stabilization  terms.  We  see 
that  this  VMS  formulation  is  a  type  of  Petrov-Galerkin  method. 

2.5  VMS  Formulation  for  Navier-Stokes 

We  apply  the  VMS  method  derived  above  to  the  Navier-Stokes  equa¬ 
tions.  This  follows  a  formulation  by  Kees  et  ah  [47].  We  assume  that  all  of 
the  test  and  trial  spaces  can  be  written  as  a  direct  sum  of  the  discrete  hnite 
element  spaces  introduced  above  and  an  inhnite  dimensional  space  of  functions 
describing  the  subgrid  scale  behavior.  So,  ©  V’',  W’'  =  ©  W’', 

=  Vl  and  and  all  of  the  functions  can  be  writ¬ 

ten  as  a  unique  sum  of  functions  in  the  grid  scale  and  subgrid  scale  spaces 
v/j  +  =  v’'  G  V’',  =  w’'  G  ,  ph  +  p  =  p  G  V^,  and 

wf^  +  vjP  =  vjP  E  W^. 

The  adjoint  operators  required  for  the  VMS  formulation  are 

L:,r><  =  (2.54) 

Ll,v^h  =  -VwW;]  -  z/Aw]]  (2.55) 


20 


41 


L!  w"  = 


dwly  dwl^ 


(2.56) 


^  \  dx  '  dy  '  dz 

where  L*  ^  and  are  the  adjoint  operators  corresponding  to  the  divergence 
and  pressure  gradient  operators.  is  obtained  by  linearizing  the  fluxes  of 
the  momentum  equation. 

Let  v/i  jj  denote  at  time  step  n.  The  ASGS  approximations  at  time 
step  n  are  given  by 

P  =  TpRp  (2.57) 

V  =  r^R,,  (2.58) 


where 


Tp  =  4z/  +  2p||v/j,„_i||/i  +  \D[\h‘^ 
1 


r„  = 


and 


4p/h2  +  +  |Dc 


-Rp  =  V  •  v/j 


=  DiVft  +  ■  Vv/j  -  z/Avft  -  g  +  -Vp, 

P 


(2.59) 

(2.60) 

(2.61) 

(2.62) 


and  h  =  he  for  k  =  1  and  h  =  hel‘1  for  k  =  2.  Notice  that  we  lag  the 
stabilization  parameters  in  time. 

The  VMS  finite  element  formulation  is:  And  ph  G  V)f  and  v/j  G  such 

that 


Yh-Vwldn+  yL*  wldn  =  -  hPdT,  'ivjPeWP  (2.63) 
In  Jn  ’  Jr?, 


21 


42 


dvh 

dt 


wr  -  Vft  (g)  Vft  ■  VwldQ  = 


/  z/(x)V(v/j  +  v^)  ■ 

'n 

[  pL*p,v^ldn 


h^^dT  ywi  e  wi 


(2.64) 


2.6  Turbulence  Modeling 

Assuming  a  porous  domain  with  the  mean  grain  size  as  the  length 
scale,  the  effects  of  turbulence  in  incompressible  flow  begin  to  be  noticeable 
at  Reynolds  numbers  near  70.  When  Re  nears  100  there  is  turbulent  flow 
in  about  half  of  the  domain.  Near  200  there  is  turbulence  everywhere  in  the 
flow,  and  above  1000  the  flow  has  fully  developed  turbulence.  A  great  deal 
of  energy  from  the  mean  flow  is  transformed  into  turbulent  kinetic  energy 
which  cascades  through  eddies  of  decreasing  size  until  it  dissipates  from  the 
system.  In  most  cases  it  is  unreasonable  to  resolve  a  computational  domain 
down  to  the  smallest  length  scales  of  turbulent  eddies,  so  it  becomes  essential 
to  incorporate  a  turbulence  model  to  capture  the  effects  of  the  subgrid  scale 
turbulent  eddies.  There  are  two  main  types  of  turbulence  models:  Reynolds 
Averaged  Navier-Stokes  (RANS)  and  Large  Eddy  Simulation  (LES).  We  focus 
on  LES  methods. 
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2.6.1  Large  Eddy  Simulation 

The  LES  equations  are  a  filtered  version  of  the  Navier-Stokes  equations 


(2.65) 


dui  ^  duiUj  1  dp  ^  d'^Ui  drij 

dt  dxi  pdxi  dxidxi  dx-i 


dt  dxj  pdxi  dxjdxj  dxj 

The  bar  signihes  a  hltered  term.  An  LES  hltering  operation  is  dehned  by 


(2.66) 


(p{x,t)  = 


0(r,  t')G{x  —  r,t  —  t')dt'dr 


(2.67) 


where  G  is  the  hlter  convolution  kernel.  The  choice  of  G  dehnes  the  type  of 
LES.  (j)  can  be  separated  into  a  hltered  component  0  and  sub-hltered  compo¬ 
nent  0',  so  0  =  0-f0'.  In  order  to  close  this  system  (2.66),  we  must  model  the 
term  r*,-.  The  most  common  model  for  this  term,  is  the  Smagorinsky  model; 


Tij  —  — 2z/'rS'j. 


(2.68) 


where 


Q  -  1  I  ^ 

2  dxi 


(2.69) 


is  an  entry  of  the  strain  rate  tensor.  The  Smagorinksy  model  is  in  a  class  of 
LES  models  called  eddy  viscosity  models.  We  can  incorporate  the  above  into 
an  eddy  viscosity  with 


ut  —  {Gs^gY  J2SijSij, 


(2.70) 


where  A„  is  the  grid  size  and  Gg  is  the  Smagorinksy  Constant  taken  as 


Gg  =  0.17. 
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The  relationship  between  the  Smagorinsky  model  for  Tij  and  the  correspond¬ 
ing  Smagorinsky  LES  hlter  is  based  on  the  energy  spectrum  of  homogenous, 
isotropic  turbulence. 


Use  of  the  Smagorinksy  LES  model  can  be  problematic  if  the  flow  is 
not  fully  turbulent  everywhere  in  the  flow.  In  this  case,  it  tends  to  be  overly 
diffusive  and  can  cause  unphysical  results.  The  Germano  dynamic  model  [31] 
is  an  LES  model  that  does  not  assume  that  Cg  is  constant.  Cg  varies  spatially 
and  temporally  depending  on  local  flow  characteristics.  This  allows  it  to  be 
used  accurately  for  flows  with  varying  amounts  of  turbulence  throughout  the 
flow.  It  involves  using  a  second  larger  hlter  "  called  the  test  hlter.  The  resolved 
stress  tensor  C  is  dehned  as 

C  =  Tr-  fl,  (2.71) 

where  =  tUuJ  —  UiUj  is  the  residual  stress  tensor  for  the  test  hlter  scale 
and  f[j  =  UiUj  —  uiUj  is  the  residual  stress  tensor  for  the  grid  hlter.  C  can  be 
considered  the  contribution  to  the  subgrid  scale  stresses  at  the  length  scales 
between  the  two  hlter  sizes.  The  resulting  equation  for  Cg  is 


(j2  _ 

"  MijMij 


(2.72) 


where  M.ij 


2A^ 


and  a 


regions  gives 


^2  _ 


=  .  Averaging  over  small 

(2.73) 
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A  full  analysis  of  LES  models  is  found  in  [80].  Now  our  system  is  the  Navier- 
Stokes  equations  with  an  added  eddy  viscosity 

(9v  Vn 

—  +  V-(v(g)v-(z/  +  z/t) Vv)  H - g  =  0.  (2.74) 

ot  p 

The  rate  of  strain  should  be  small  for  low  Reynolds  number  flows,  so  vt  will 
only  have  a  signihcant  effect  for  higher  Reynolds  nnmber  flows. 

Yoke  [107]  and  Menevean  and  Lnnd  [69]  developed  more  practical  meth¬ 
ods  for  calculating  Cg-  By  assnming  a  Pao  tnrbnlent  energy  spectrum,  for  fully 
laminar  flows  we  take  Cg  =  0,  and  for  transitional  and  turbulent  flows,  we  take 

Cl{Re^)  =  0.027  x  (2.75) 

where  Rca  =  jg  the  mesh-Reynolds  nnmber  at  a  given  point  in  the 

flow.  This  approximation  to  the  fnll  dynamic  Smagorinsky  method  is  the  one 
that  we  ntilize. 

2.6.2  LES  Verification 

In  order  to  establish  that  the  LES  model  presented  above  is  a  realistic 
one  we  evalnate  it  against  experimental  data.  Lin  et  al.  [59]  performed  an 
experimental  study  of  flow  through  rigid  vegetation.  They  simnlated  a  vege¬ 
tated  held  with  several  different  periodic  arrangements  of  rigid  dowels.  They 
experimented  with  both  emergent  and  snbmerged  vegetation  and  smooth  and 
rongh  beds.  We  validated  the  LES  method  by  comparing  to  their  Experiment 
1.2.  In  this  setnp,  rigid,  smooth  dowels  0.076  m  tall  and  0.00635  m  in  diameter 
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are  arranged  in  a  grid  with  their  centers  0.1016  m  apart.  Their  flume  was  4.3 
m  long  and  0.3  m  wide  with  a  constant  slope  of  0.003  and  the  dowels  started 
1.3  m  downstream.  The  flow  rate  was  0.0057  m^/s  and  the  flow  depth  was 
0.071  m.  Therefore  the  dowel  is  emergent.  The  flow  is  fully  developed  by  the 
time  it  is  2.25  m  downstream  of  the  inflow  boundary.  They  made  instanta¬ 
neous  longitudinal  velocity  measurements  with  one-dimensional  laser  Doppler 
velocimeters  at  specific  locations  in  the  flow  and  throughout  the  vertical  pro¬ 
file  at  these  points.  We  focus  on  3  of  these  locations.  Measurement  1  is  taken 
4  dowel  diameters  directly  downstream  from  the  center  of  a  dowel.  Measure¬ 
ment  2  is  taken  8  dowel  diameters  directly  downstream  from  the  center  of  the 
dowel.  Measurement  4  is  taken  8  dowel  diameter  directly  downstream  and  8 
dowel  diameters  to  the  right  (with  respect  to  flow  direction)  of  the  center  of 
the  dowel. 

It  would  be  unrealistic  for  us  to  fully  simulate  their  complete  setup.  We 
simulate  on  a  0.0762  m  x  0.0762  m  x  0.071  m  rectangular  prism  with  a  cylinder 
of  diameter  0.00635  mm  at  (0.0254  m,  0.0254  m).  Our  mesh  has  223,930 
unstructured  tetrahedral  elements.  We  use  standard  parameters  p  =  998.0 
kg/m^  and  u  =  1.0  x  10“®  m^/s  for  water.  We  impose  no  slip  boundary 
conditions  on  the  bottom  and  on  the  cylinder,  a  free  surface  condition  on 
the  top,  periodic  conditions  on  the  transverse  sides,  outflow  conditions  on  the 
back,  and  a  Dirichlet  condition  on  the  front  face.  Starting  with  zero  initial 
conditions  we  ramp  up  the  velocity  over  time  up  to  0.2676  m/s,  which  is  the 
mean  velocity  of  the  flow  in  the  experiment.  We  then  let  the  flow  develop 
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over  10  s  and  take  instantaneous  snapshots  of  the  velocities  prohles  at  the  3 
measurement  locations  shown  in  Figure  2.1.  These  qualitatively  match  quite 
well  with  Figure  2  in  the  paper  of  Liu  et  al.  [59].  As  with  their  results, 
the  general  trend  is  that  the  lowest  velocities  occur  at  locations  behind  the 
cylinder.  We  also  see  that  the  velocity  at  location  4,  in  the  free  stream,  is 
always  the  highest.  As  with  the  experiments,  we  see  that  the  dowel  has  a 
large  effect  on  the  prohle,  with  the  velocity  prohles  being  nearly  vertical  in 
the  intermediate  depths.  As  with  their  results,  we  see  spikes  in  velocity  just 
above  the  bed.  This  is  due  to  a  horseshoe  vortex  at  the  base  of  the  cylinder. 
Capturing  this  spike  validates  our  method.  The  magnitude  of  the  velocities 
is  very  close  to  the  experimental  results  at  all  3  locations.  This  shows  that 
the  LES  method  can  realistically  model  the  complex  how  dynamics  for  high- 
velocity  how  around  these  types  of  structures. 
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Measurement  Location  1 


Figure  2.1:  Longitudinal  velocity  profiles  with  respect  to  depeth  (z)  at  loca¬ 
tions  1,  2  and  4  respectively  for  LES  verification  study. 
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Chapter  3 

Computational  Upscaling  of  Flow  Through 
Domains  Containing  Rigid  Vegetation 

3.1  Introduction 

For  many  types  of  physical  processes  there  are  different  characteristics 
at  different  length  scales.  It  is  often  computationally  unreasonable  to  resolve 
the  smallest  length  scales  of  interest  for  a  large-scale  problem.  Upscaling  is  the 
process  of  analyzing  small-scale  dynamics  to  describe  larger  scale  phenomena. 
Flow  through  vegetated  regions  is  an  example  of  one  of  these  processes  where 
the  dynamics  are  at  many  length  scales.  For  a  large  vegetated  region,  it 
would  be  nearly  impossible  to  resolve  all  of  the  complex  geometry  and  flow 
characteristics,  so  we  must  use  upscaling. 

Mathematical  homogenization  is  a  common  upscaling  technique.  It  in¬ 
volves  mathematically  deriving  effective  equations  to  describe  large-scale  dy¬ 
namics  by  satisfying  to  some  degree  the  dynamics  of  the  smaller  scales.  Most 
homogenization  theory  exists  for  linear  problems.  The  topic  of  homogeniza¬ 
tion  for  nonlinear  problems  is  mostly  open,  and  most  work  in  the  held  involves 
problems  that  are  only  slightly  nonlinear.  We  use  techniques  from  mathemat¬ 
ical  homogenization  for  highly  viscous  hows;  however,  for  inertial  hows,  which 


29 


50 


Figure  3.1:  The  periodic  cell  Y . 

we  are  more  likely  to  see  in  nature  in  vegetated  domains,  we  must  perform 
computational  upscaling.  This  chapter  explains  a  method  of  computational 
upscaling  that  we  have  published  in  [65]. 

3.2  The  Concept  of  Drag 

Fluid  viscosity  causes  resistance  in  the  flow  around  immersed  bodies. 
Viscous  effects  can  produce  three  different  types  of  resistance  as  described  by 
Rouse  [90].  At  low  Reynolds  numbers,  inertial  effects  of  flow  are  negligible 
compared  to  those  caused  by  viscous  stress.  We  see  that  these  viscous  effects 
extend  a  great  distance  into  the  surrounding  flow,  causing  a  widespread  dis¬ 
tortion  of  the  flow  pattern.  This  is  known  as  “deformation  drag.”  At  higher 
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Reynolds  numbers,  the  deformed  region  of  the  flow  is  much  smaller,  limited  to 
a  thin  layer  surrounding  the  body.  Therefore  the  resulting  shear  is  along  the 
boundary  surface.  We  call  this  “surface  drag.”  If  the  form  of  the  body  is  such 
that  there  is  separation  in  the  flow,  it  produces  “form  drag.”  In  this  scenario, 
there  is  a  lower  intensity  of  pressure  in  the  wake  of  the  object  which  causes  a 
resultant  force  which  opposes  the  motion.  Under  certain  conditions  form  drag 
can  reduce  viscous  stresses  to  insignihcant  values. 

We  write  the  relationship  for  the  force  of  drag  opposing  motion  as 

F  =  C,A^  (3.1) 


where  F  is  the  drag  force,  A  is  the  cross  sectional  area  of  the  body,  V  is  the 
magnitude  of  the  velocity  of  the  flow,  and  Cd  is  a  drag  coefficient.  We  note 


that 


Cd  = 


(3,2) 


ApV‘^12 

is  dependent  on  Reynolds  number.  This  coefficient  is  similar  to  the  Darcy 
friction  factor  for  the  Darcy- Weisbach  equation  used  in  modeling  flow  in  rough 
pipes.  In  multiple  dimensions,  we  can  write  the  drag  equation  as 


F 


CdA 


p|V|V 

2 


(3.3) 


where  now  Cd  is  a  tensor  and  V  is  the  mean  velocity  vector.  For  calculating 
Cd,  we  take  F  as  the  hydraulic  gradient 


F  =  -LAiVp  -  pf) 


(3.4) 
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where  L  is  the  domain  length  and  f  is  a  forcing  term  which  includes  the 
gravitational  constant.  In  this  paper,  we  study  how  Cd  changes  over  a  large 
range  of  Reynolds  numbers  for  various  vegetated  domains. 


3.3  Upscaling  Laws 


We  will  use  some  techniques  from  the  study  of  flow  through  porous 
media  for  our  analysis  of  flow  through  vegetated  regions.  Tightly  packed 
vegetation  can  be  analogous  to  a  porous  medium.  We  dehne  Reynolds  number 
as 


V 


(3.5) 


where  V  is  the  mean  flow  rate  over  the  total  volume  including  the  obstructions 
(specihc  discharge),  d  is  the  mean  diameter  of  the  individual  plants,  and  z/  is 
the  kinematic  viscosity  of  the  fluid.  Reynolds  number  is  a  non-dimensional 
value  representing  the  ratio  of  inertial  effects  to  viscous  effects.  In  most  studies 
of  subsurface  flow,  it  is  assumed  that  viscous  effects  dominate,  so  inertial  effects 
can  be  ignored. 


3.3.1  Darcy’s  Law 

For  very  low  Reynolds  number  flows  {Re  <  1)  Darcy’s  law  is  the  basic 
constitutive  equation  for  flow  in  porous  media.  It  was  originally  determined 
experimentally  by  Henry  Darcy  in  1856.  As  dehned  by  Bear  [5],  Darcy’s 
law  says  that  for  a  homogeneous  incompressible  fluid,  the  three-dimensional 
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constitutive  equation  is 

v  =  -  — (Vp-pf),  (3.6) 

p 

where  v  is  the  mean  velocity  vector,  Vp  is  the  pressure  gradient,  K  is  the 
“hydraulic  permeability”  tensor  with  units  p  is  the  dynamic  viscosity  of 
the  fluid,  p  is  the  fluid  density,  and  f  is  the  forcing  vector  containing  the  grav¬ 
itational  constant.  The  hydraulic  permeability  tensor  depends  solely  on  the 
properties  of  the  porous  medium.  Since  we  are  dealing  with  an  incompressible 
fluid,  we  also  assume  that  the  flow  is  divergence  free: 

V-v  =  0.  (3.7) 

For  Reynolds  numbers  small  enough  for  Darcy’s  law  to  be  valid,  we 
only  need  to  know  the  hydraulic  permeability  tensor  to  be  able  to  accurately 
describe  the  mean  flow.  K  is  calculated  by  mathematical  homogenization. 

Assume  the  domain  is  periodic  and  perforated.  This  can  be  seen  in 
Figure  3.1.  Let  F  be  a  standard  periodic  cell  with  standard  obstacles,  G  CC  F, 
with  piecewise  smooth  boundaries  Fob^.  The  remainder  is  D  =  F  \  G.  This  is 
the  “fluid  region” .  We  call  the  outer  fluid  boundary  of  the  cell  that  does  not 
touch  an  obstacle  Tauter-  We  assume  the  standard  geometry  of  F  is  repeated 
all  throughout  R^. 

A  hxed  domain  A  is  given  by  intersecting  the  e  multiple  of  this  periodic 
geometry.  This  is  denoted  by 


D"  =  A  n  (eD) 
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and 

ri,,  =  An{erobs). 

For  the  simple  theory,  we  assume  that 

dA  n  eTobs  =  0; 

however,  in  practice  this  condition  is  not  necessary. 

Now,  we  perform  homogenization  on  the  standard  Stokes  problem  to 

formally  derive  Darcy’s  Law  as  shown  by  Hornung  [37].  We  start  with  the 

following  problem  on  the  pore  scale 

e^/iAv'^(a;)  =  Vp'^{x),  x  G 

V  •  v'^(a;)  =  0  a;  G  (3.8) 

v^(a;)  =  0,  a;  G  F^^, 

where  fi  is  the  dynamic  viscosity,  v  is  velocity  and  p  is  pressure.  We  assume 
that  the  unknowns  p'^  and  can  be  written  as  asymptotic  expansions  with 
respect  to  e.  That  is 

v"(a;)  =  vo(a;,|/)  +  evi(a;,|/)  +  e^V2(a;,  y)...  (3.9) 

and 

P%x)  =  po{x,  y)  +  epi(a;,  y)  +  e^P2{,x,  y)...  (3.10) 

where  the  coefficient  functions  Vi  and  Pi  are  Y-periodic  with  respect  to  y.  This 
results  in  the  equations 

e°p,AyVo{x,y)  +  =  e~^VyPoix,y) 

+  e^(yyPi{x,y)  +  V^po{x,y))  +  e\...)  +  ..., 
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e  ■  vo(a;,|/)  +  e°(Vj^  ■  vi(a;,|/)  +  Va:  ■  vo(a;, I/)  +  e^(.. .)  +  •••=  0,  (3.11) 

for  all  y  E  fl,  and 

e°vo(a;,  y)  +  eVi(a;,  y)  +  ...  =  0,  (3.12) 

for  all  y  eT.  Now  we  compare  coefficients.  The  e~^  term  in  the  first  equation 
gives 

'^yPo{x,y)  =  0,  (3.13) 

for  all  y  E  fl,  hence,  po{x,y)  =  po{x)  independently  of  y.  The  e°  term  yields 

pAyVoix,  y)  =  WyPi{x,  y)  +  (3.14) 

for  dW  y  eVL  and  the  term  in  the  second  equation  gives 

V  ■  vo(a;,|/)  =  0,  (3.15) 

for  all  y  G  fl.  As  usual 

VxPo(a:)  =  '^ejdjPo{x)  (3.16) 

j 

where  ej  is  the  unit  vector  in  the  j  direction. 

Now  the  cell  problem  is  to  hud  the  Y-periodic  vector  helds  Wj(|/)  with 
components  Wji{y)  that  solve  the  Stokes  problems, 

-^j  y 

V  ■  Wj{y)  =  0  yEVL  (3-17) 

Wj(|/)  =  0  yETobs 

where  'Kj{y)  are  the  corresponding  Y-periodic  pressure  helds,  Gj  is  the  unit 
vector  in  the  j  direction,  v  is  the  kinematic  viscosity  and  p  is  the  huid  density. 
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Notice  that  the  third  equation  is  our  standard  no-slip  condition  on  the  walls 
of  the  vegetation.  Using  these  cell  functions,  we  easily  obtain 

Vo  (a:,  y)  =  --Y1  (3-18) 

^  j 

The  average  vector  held  is  dehned  by 


u(x)  =  /  \o{x,y)dy, 


and  its  i*^  component  is  expressed  by 


Ui{x)  = - ^  kijd^.po{x), 

J 

with 

kij  =  /  Wji{y)dy.  (3.19) 

Jn 

By  introducing  the  tensor  K  =  k^j  we  get  the  short  expression 

_  1 

u(a;)  = - KVpo(a^) 

/i 

which  is  Darcy’s  Law. 

We  must  still  show  that  u  is  divergence  free.  From  the  e°  term  of  (3.11). 


Vy  ■  vi{x,y)  +  Vx  ■  vo{x,y)  =  0, 
for  al  y  G  D.  By  integrating  over  D, 
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v^-ui^x)  =  /  V-\o{x,y)dy 

Jn 

=  -  [  Vy-vi{x,y)dy 


=  -  n-vi{x,y)dTouter{y) 

r outer 

=  -  n-Vi{x,y)dTouter{y)  -  n-Vi{x,y)dTouter{y)- 

^ outer 

The  boundary  integral  over  T outer  is  zero  due  to  the  term  with  in  (3.12) 
and  the  boundary  integral  over  dY  is  zero  due  to  the  Y-periodicity  of  vi(a;,  y). 
Therefore  u  is  divergence  free. 

We  can  also  write  Darcy’s  Law  in  terms  of  the  multi-dimensional  drag 
coefficient.  Using  Darcy’s  Law  and  (3.4)  we  have 


F  =  LAyK^^V. 

Comparing  with  the  drag  equation  we  obtain 

2Ly 


Cd  = 


p|v| 


K\ 


(3.20) 


3.3.2  Darcy-Forchheimer  Flow 

It  is  observed  experimentally  that  for  flows  with  Reynolds  numbers 
greater  than  1  the  ratio  between  hydraulic  gradient  and  velocity  gradually 
increases  with  increasing  flow  rate.  It  is  generally  thought  that  this  increase 
is  because  inertial  forces,  with  increasing  Reynolds  number,  become  much 
larger  than  the  viscous  forces.  In  this  range,  nonlinear  behavior  is  seen.  The 
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inertial  terms  in  Navier-Stokes  are  no  longer  negligible,  so  the  homogenization 
of  Stokes  equation  is  no  longer  the  correct  upscaling.  Forchheimer  noticed  from 
empirical  results  that  for  higher  Reynolds  numbers,  the  relationship  resembles 
Darcy’s  Law  plus  a  quadratic  term.  This  is  called  the  Darcy-Forchheimer 
equation 

-(Vp-pf)  =/iK”^V  +  /3p|V|V.  (3.21) 

The  non-Darcy  coefficient  (3  is  usually  taken  as  a  constant,  but  it  is  likely 
a  tensor  for  anisotropic  media.  It  is  generally  calculated  experimentally  by 
htting  to  data.  While  the  Darcy-Forchheimer  equation  has  matched  data  from 
a  large  number  of  experiments  quite  well,  it  is  often  unacceptable,  especially 
over  a  large  range  of  Reynolds  numbers  as  explained  in  Bear  [5]. 

3.3.3  Other  non-Darcy  models 

There  have  been  many  attempts  at  deriving  the  Darcy-Forchheimer 
Law  using  mathematical  homogenization.  One  of  the  hrst  was  by  Lions  [58]. 
Similar  techniques  were  done  by  Sanchez  [91].  These  types  of  techniques  found 
that  the  obtained  homogenization  problems  lead  to  polynomial  filtration  laws. 
This  has  been  studied  by  Mei  and  Auriault  [67],  Wodie  and  Levy  [112],  and 
Rasoloarijaona  and  Ariault  [83].  The  results  of  these  studies  sees  quadratic 
terms  canceling,  giving  a  cubic  law 


-(Vp-pf)  =/iK-^V  +  /3p|V|2v.  (3.22) 
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This  has  been  verihed  numerically  for  very  small  Reynolds  numbers  by  solving 
the  Navier-Stokes  equations  by  Coulaud  et  ah  [14],  Rojas  and  Koplik  [89],  and 
Koch  and  Ladd  [50].  However,  this  has  only  been  verihed  where  Darcy’s  Law 
is  approximately  valid,  and  the  quadratic  law  seems  to  be  more  appropriate  for 
higher  Reynolds  number  hows.  Some  methods  provide  a  result  with  quadratic 
and  cubic  terms 

-(Vp  -  pf)  =  +  /3p|V|V  +  ypiVpV.  (3.23) 

Some,  like  Scheidegger  [93],  have  also  used  power  laws  of  the  form, 

-(Vp  -  pf)  =  pK-^V  +  (3.24) 

to  ht  empirical  data.  For  (3.24)  and  (3.23),  /?,  7,  ai,  and  02  are  model  pa¬ 
rameters  that  must  be  found  empirically. 

3.3.4  Higher  Reynolds  Number  Flows 

From  experimental  data  we  know  that  as  Reynolds  number  increases 
above  the  region  where  Darcy’s  Law  is  valid,  the  nature  of  the  how  changes. 
At  low  Reynolds  numbers  we  see  deformation  drag;  streamlines  are  greatly 
ahected  by  the  obstacles.  As  Reynolds  number  increases,  we  see  separation 
occur  and  form  drag  occurs  with  less  deformation.  Eventually,  streamlines 
shift  and  hxed  eddies  form  in  the  wake  of  obstacles.  The  size  of  the  eddies 
increases  as  Re  increases.  Around  Re  =  70  turbulence  begins  to  occur.  When 
Re  reaches  75  the  how  is  turbulent  in  approximately  half  of  the  domain,  and 
when  Re  nears  200  there  is  turbulence  everywhere  in  the  how.  This  analysis 
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is  discussed  at  length  by  Bear  [5].  At  high  Reynolds  numbers  we  see  the  drag 
force  begin  to  balance  with  the  driving  force.  Therefore,  we  see  Cd  become 
constant.  Chaudhary  et  al.  [12]  explain  that  the  change  in  drag  from  low 
Reynolds  numbers  to  moderate  Reynolds  numbers  is  due  to  the  increasing 
size  of  eddies  within  the  pores.  For  high  Reynolds  numbers,  the  eddies  £11  up 
almost  all  of  the  pore  space  and  can  no  longer  grow,  so  hydraulic  conductivity 
and  drag  coefficient  become  nearly  constant. 

Darcy’s  Law  and  the  derivation  of  hydraulic  conductivity  above  do  not 
accurately  model  the  nonlinear  behavior  in  this  range.  As  we  see  in  the  previ¬ 
ous  sections,  the  Darcy-Forchheimer  Law  and  the  related  polynomial  formulas 
do  have  some  theoretical  arguments  and  seem  to  match  experimental  data  for 
the  initial  nonlinear  region.  However,  there  is  not  a  theoretical  basis  to  ex¬ 
trapolate  these  results  to  higher  Reynolds  numbers.  In  subsurface  modeling, 
it  is  uncommon  to  deal  with  Reynolds  numbers  this  large.  Our  main  goal  is 
to  analyze  the  drag  in  this  transitional  region  between  Stokes  and  turbulent 
flow  so  that  it  can  be  applied  to  flow  through  vegetated  regions  where  we  do 
see  Reynolds  numbers  in  this  range. 

3.4  Method 

3.4.1  Calculation  of  Hydraulic  Conductivities 

In  order  to  describe  Cd  for  low  Reynolds  number  flows  we  can  use  the 
above  method  of  mathematical  homogenization  to  find  hydraulic  conductiv¬ 
ities  for  various  2D  and  3D  domains.  We  computationally  solve  the  Stokes 
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Figure  3.2:  3D  domains  containing  cylinders  and  cones. 

problems,  (3.17).  Since  we  are  solving  on  a  periodic  cell,  we  use  periodic 
boundary  conditions.  The  weak  formulation  of  these  problems  is  a  simpler 
version  of  what  is  done  in  Chapter  2.  Once  the  periodic  w  and  vr  helds  are 
found,  we  use  (3.19)  to  calculate  K,  integrating  with  numerical  quadrature. 
Then  using  (3.20),  we  are  able  to  calculate  the  drag  coefficients  for  these  vis¬ 
cosity  dominated  flows. 

For  higher  Reynolds  number  flows,  we  must  perform  “computational 
upscaling.”  Instead  of  solving  the  linear  Stokes  problems  (3.17),  we  solve  the 
cell  LES  problem: 

^  V  ■  (wj  0  wj  -  (z/  z/r)  Vwj)  -H^-fj=0  y  eVt 

V  ■  Wj(|/)  =  0  1/  G  D  (3.25) 
Wj(2/)  =  0  yeVohs 

for  a  large  range  of  on  the  periodic  cell  to  “steady-state”.  For  moderate 
Reynolds  numbers  {Re  <  100),  there  is  a  true  steady-state  solution  that  can  be 
found  using  the  steady  Navier-Stokes  equations  which  ignore  the  acceleration 
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term.  For  higher  Reynolds  numbers,  the  full  cell  LES  problem  (3.25)  must 
be  solved  out  in  time  until  the  solution  becomes  essentially  steady.  When  the 
change  of  the  volume  averaged  velocity  is  within  a  small  tolerance  between 
multiple  time  steps,  we  consider  the  mean  flow  to  be  essentially  steady.  For 
very  high  Reynolds  numbers,  the  flow  never  becomes  steady,  so  we  must  time- 
average  the  solution. 

After  solving  the  cell  LES  problem,  we  calculate  corresponding  drag 
forces  and  volume-averaged  velocities.  As  for  extremely  low  Reynolds  number 
flows,  we  solve  the  Stokes  cell  problems  from  mathematical  homogenization  to 
calculate  K.  We  perform  least-squares  hts  to  calculate  the  non-Darcy  param¬ 
eters  that  best  ht  the  data  over  the  whole  range  of  Reynolds  numbers.  Then, 
we  can  parametrize  Cd  using  Darcy  and  non-Darcy  laws. 

3.4.2  Domains  and  Boundary  Conditions 

For  simplicity,  for  2D  computations  we  use  a  1  m  x  1  m  box  for  the 
porous  periodic  cell  as  seen  in  Figure  3.1.  For  3D,  we  utilize  2  types  of  do¬ 
mains  that  are  both  Imxlmxlm  cubes  as  seen  in  Figure  3.2.  For  our 
3D  simulations  we  assume  periodicity  in  the  2  horizontal  directions.  The  hrst 
3D  domain  contains  a  cylindrical  obstacle  penetrating  the  entire  domain.  This 
will  be  used  to  simulate  “emergent  vegetation”:  vegetation  that  rises  above 
the  water  level.  The  second  3D  domain  seen  in  Figure  3.2  contains  a  conical 
obstacle  which  does  not  reach  the  top  of  the  domain.  This  models  “submerged 
vegetation,”  that  which  lies  completely  below  the  water  level.  We  call  the  re- 
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Simulation  type 

Obstacle  Diameter(m) 

Number  of  Elements 

2D  Porous  Media 

0.5 

4014 

2D  Porous  Media 

0.55 

4092 

2D  Porous  Media 

0.6 

3934 

2D  Porous  Media 

0.65 

3918 

3D  Emergent  Cylinder 

0.6 

259,151 

3D  Emergent  Cylinder 

0.7 

219,923 

3D  Submerged  Cone 

0.35 

339,432 

Table  3.1;  Number  of  elements  for  all  simulations. 


gion  of  the  cell  containing  the  fluid  hi.  The  boundary  of  hi,  T,  can  be  separated 
as  T  =  TperjorficUrtopUrboitomUTofe^  where  Tperiodic  contains  the  front,  back,  left, 
and  right  faces  of  the  cubic  domain.  No-slip  boundary  conditions  are  imposed 
on  Tobs  U  r bottom,  periodic  boundary  conditions  are  imposed  on  Tperiodic  and  a 
free  surface  boundary  condition  is  imposed  on  T top- 

3.4.3  Numerical  Method 

We  use  a  locally  conservative,  stabilized  hnite  element  method,  outlined 
in  [47]  and  Chapter  2,  to  hnd  weak  solutions  to  (3.17)  and  (3.25)  with  zero 
initial  conditions  for  pressure  and  velocity.  For  higher  Re  cases  we  use  implicit 
time  stepping  until  the  mean  velocities  reach  a  stationary  state.  We  utilize 
linear  Lagrangian  basis  functions  on  unstructured  triangular  and  tetrahedral 
meshes.  Table  3.1  gives  the  number  of  elements  for  each  simulation.  When  we 
include  the  time  derivative,  we  convert  the  initial-boundary  value  problem  into 
a  sequence  of  boundary  value  problems.  All  simulations  were  run  in  parallel 
on  Texas  Advanced  Computing  Center  (TACC)  machines. 
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(a)  Re  =  0.09 


(b)  Re  =  9.4 


(c)  Re  =  34.9 


Figure  3.3:  Streamlines  (left)  and  pressure  contours  (right)  for  instantaneous 
snapshots  of  flows  of  various  Reynolds  numbers  in  2D  porous  media  domain 
with  circle  diameter  0.25  m.  Streamline  plots  have  background  coloring  show¬ 
ing  velocity  magnitude  and  pressure  contour  plots  have  background  coloring 
showing  pressure.  Red  coloring  means  higher  values  and  blue  means  smaller. 
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(a)  i?e  =  0.1 


(b)  Re  =  9.8 


(d)  Re  =  253.1 


(e)  Re  =  1128.1 


r  % 

X 

X 

X 


Figure  3.4:  Streamlines  (left)  and  pressure  contours  (right)  for  instantaneous 
snapshots  of  flows  of  various  Reynolds  numbers  in  2D  porous  media  domain 
with  circle  diameter  0.6  m.  Streamline  plots  have  background  coloring  show¬ 
ing  velocity  magnitude  and  pressure  contour  plots  have  background  coloring 
showing  pressure.  Red  coloring  means  higher  values  and  blue  means  smaller. 
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3.5  Results 


Our  analyses  involve  the  change  in  behavior  of  flow  over  a  range  of 
Reynolds  numbers.  We  take  as  our  “true”  values  the  results  of  solving  (3.25). 
For  simplicity,  we  assume  p  =  1  kg/m^  and  v  =  1  m^/s.  This  is  acceptable 
because  we  use  the  non-dimensional  Reynolds  number  as  our  main  variable 
representing  flow  speed.  The  hlter  size  is  taken  as  the  diameter  of  the  cell. 
The  flows  are  driven  by  the  forcing  term  f  in  (3.25).  We  then  recover  velocity 
fields  and  pressure  values. 

3.5.1  Qualitative  Analysis  of  Drag 

The  nature  of  the  flow  and  drag  change  as  Reynolds  number  increases. 
We  solve  (3.25)  driven  by  a  range  of  forcing  values  f.  In  Figure  3.3  and  Figure 
3.4  we  show  streamlines  and  pressure  contours  for  2D  flows  with  obstacles 
of  diameter  0.25  m  and  0.6  m  respectively.  These  domains  are  our  standard 
porous  media  domains  of  arrays  of  circles. 

In  Figure  3.3  (a)  with  Re  =  0.09  we  see  Stokes  flow.  The  flow  is  almost 
completely  dominated  by  viscous  forces.  We  see  the  viscous  effects  extend¬ 
ing  far  into  the  surrounding  fluid.  Notice  that  the  streamlines  are  deformed 
throughout  the  entire  flow.  This  is  the  “deformation  drag”  discussed  in  Sec¬ 
tion  2.1.  The  pressure  held  is  completely  symmetric.  The  increase  in  pressure 
on  the  front  of  the  obstacles  is  equal  and  in  the  same  pattern  as  the  decrease 
in  pressure  on  the  back  side.  We  expect  this  for  Stokes  how.  In  (b)  with 
Re  =  9.4,  we  see  that  while  there  is  still  a  fair  amount  of  deformation  to  the 


46 


67 


flow,  there  is  much  less  than  in  (a).  This  signihes  that  the  flow  is  no  longer 
dominated  by  viscous  forces.  There  are  now  nonlinear  inertial  effects.  Since 
flow  separation  is  not  occurring  here  this  is  considered  to  be  “surface  drag.” 
The  pressure  contour  patterns  are  no  longer  symmetric.  In  (c)  with  Re  =  34.9, 
we  see  that  laminar  recirculation  has  started  behind  the  obstacles.  Note  that 
they  are  not  turbulent  eddies,  as  turbulence  does  not  occur  at  this  Reynolds 
number.  There  is  very  little  deformation  to  the  drag  in  the  region  between  the 
obstacles.  The  drag  that  we  see  here  is  “form  drag.”  In  (d),  with  Re  =  110.0 
we  see  that  the  recirculation  regions  have  grown  greatly  in  size,  as  we  expect. 
There  is  even  less  flow  deformation  due  to  the  obstacles.  This  flow  has  some 
turbulence,  but  not  enough  to  have  significant  effects  on  the  flow  pattern.  In 
(e)  with  Re  =  1082.5,  we  see  fully  turbulent  flow.  There  is  no  steady-state 
flow  pattern.  The  obstacles  cause  almost  no  deformation  in  the  flow  pattern. 

In  Figure  3.4  we  see  results  for  a  domain  with  circles  of  diameter  0.6  m. 
Qualitatively,  these  results  are  quite  similar  to  those  for  the  previous  domain. 
Because  of  the  much  smaller  porosity;  however,  we  see  that  the  deformation 
in  the  flow  remains  for  higher  Re.  Also  notice  the  difference  in  the  pattern 
of  the  pressure  contours  as  compared  to  the  previous  example.  In  the  more 
densely  packed  domain,  the  pressnre  contours  go  between  separate  obstacles, 
and  they  do  not  in  the  less  tightly  packed  domain.  This  resnlt  is  similar  to 
the  resnlts  from  Koch  and  Ladd  [50]  and  demonstrates  that  the  tightness  of 
particle  packing  has  a  snbstantial  effect  on  important  flow  characteristics. 
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Circle  Diameter(m) 

Separation(m) 

K(m") 

0.65 

0.35 

0.00008548 

0.6 

0.4 

0.0004262 

0.55 

0.45 

0.001148 

0.5 

0.5 

0.002362 

0.45 

0.55 

0.004175 

0.4 

0.6 

0.006701 

0.35 

0.65 

0.01007 

0.3 

0.7 

0.01449 

0.25 

0.75 

0.02024 

Table  3.2:  Values  for  hydraulic  conductivity  K  for  Porous  Media  Cells  of 
Various  Diameters 


3.5.2  Homogenization  and  Darcy’s  Law 

We  calculated  hydraulic  conductivities  K  by  solving  (3.17)  on  various 
domains.  We  use  a  highly  rehned  mesh  with  h  =  1/120  m  where  h  is  the 
maximum  element  edge  length.  For  our  standard  porous  media  domain(as 
seen  in  Figures  3.3  and  3.4)  with  varying  diameters  of  circles  we  show  the 
values  of  K  in  Table  3.2.  Graphically,  we  show  this  in  Figure  3.5. 

In  Figure  3.6  we  see,  as  we  expect,  that  Darcy’s  Law  using  the  K 
calculated  from  homogenization  is  valid  for  small  Reynolds  numbers  (<  0.25). 
Above  this  range,  we  see  inertial  effects  having  a  major  effect  on  the  flow. 
Obviously,  Darcy’s  Law  is  an  unreasonable  assumption  for  higher  velocity 
flow,  the  type  we  might  observe  in  flow  through  vegetation.  We  therefore 
investigate  the  non-Darcy  models  mentioned  above. 
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Hvdraulic  Permeabilities  for  Porous  Media  Domains  with  Obstacles  of  Different  Diameters 
0.025 1 - ^ ^ ^ ^ ^ ^ ^ - 1 


0.02 


0.015 


0.01 


0.005 


S.25  0.3  0.35  0.4  0.45  0.5  0.55  0.6  0.65 

Circle  Diameter  (m) 


Figure  3.5:  Values  for  hydraulic  conductivity  K  for  porous  media  cells  of 
various  diameters. 


Homogenization  vs.  Navier-Stokes  Data 


Figure  3.6:  A  comparison  of  Darcy’s  law  with  computationally  calculated  K 
and  Navier-Stokes  data. 
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3.5.3  Analysis  of  Non-Darcy  Laws 

We  compare  our  Navier-Stokes  data  with  the  Darcy  and  non-Darcy  laws 
introduced  earlier.  After  performing  Navier-Stokes  simulations  over  a  range  of 
forcing  terms  f,  we  have  a  collection  of  data  from  which  to  £t  parameters  for 
the  non-Darcy  models  mentioned  previously.  These  fitted  parameters  change 
drastically  depending  on  the  range  of  Reynolds  numbers  over  which  the  fit  is 
performed.  The  following  calculations  were  performed  on  our  standard  porous 
media  2D  domain  with  circles  of  diameter  0.6  m. 

Figure  3.7  (a)  shows  how  (3  from  the  Darcy-Forchheimer  equation  calcu¬ 
lated  from  least  squares  fitting  changes  depending  on  the  range  over  which  the 
fit  is  done.  Ideally,  we  could  perform  this  fit  for  low  Reynolds  number  flow  and 
it  would  remain  valid  for  higher  Reynolds  number  flow;  however,  these  results 
show  that  (3  changes  as  the  range  of  Reynolds  numbers  for  fitting  is  increased. 
Figure  3.7(b)  shows  the  Forchheimer  parameter  (3  over  a  range  of  values  from 
Re  =  1,  Re  =  10,  and  Re  =  100  compared  to  the  Navier-Stokes  data.  The 
two  fits  over  the  low  ranges  diverge  significantly  from  the  Navier-Stokes  data 
at  high  Re. 

The  results  from  fitting  the  cubic  non-Darcy  eqnation  are  in  Figure 
3.8.  The  coefficient  of  the  cnbic  term  varies  greatly  depending  on  the  range 
of  Reynolds  nnmbers  over  which  the  fitting  is  performed.  Notice  that  these 
plots  vary  so  greatly  over  the  given  range  of  Reynolds  nnmbers  that  the  y-axis 
must  be  shown  on  a  logarithmic  scale.  This  shows  that  for  higher  Re  the  cubic 
law  is  not  valid.  For  low  Reynolds  nnmbers,  it  does  fit  the  data  rather  well. 
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validating  some  of  the  homogenization  theory  by  Mei  and  Auriault  [67]  and 
Rasoloarijaona  and  Auriault  [83]. 

The  power  law  non-Darcy  equation  £t  results  are  shown  in  Figure  3.9; 
(a)  shows  how  the  coefficient  13a  changes  depending  on  the  range  of  Re  over 
which  the  fit  is  performed;  (b)  similarly  shows  this  for  the  power  a.  Both  of 
these  plots  flatten  out  at  higher  Re,  with  f3a  ~  145  and  a  ~  1.87.  The  results 
from  the  lower  Re  hts  do  a  reasonable  job  at  matching  Navier-Stokes  data. 
Figure  3.9  (c)  shows  that  the  parameter  fits  up  to  Re  =  10  and  Re  =  100 
match  extremely  well  with  the  Navier-Stokes  data. 

These  non-Darcy  models  historically  are  based  on  empirical  results  with 
data  from  relatively  low  Reynolds  number  flows  {Re  <  10).  Our  results  show 
that  they  all  do  a  reasonable  job  at  matching  the  Navier-Stokes  data  in  this  low 
Re  range;  however,  parameters  estimated  on  low  Re  data  do  not  sufficiently 
model  higher  Re  behavior.  From  our  fits  we  see  that  if  the  parameters  are 
calculated  using  a  larger  range  of  Reynolds  number,  that  they  adequately  fit 
for  higher  Re  behavior  as  well  as  lower. 

3.5.4  Cd  and  Higher  Reynolds  Number  Flows 

The  Navier-Stokes  data  can  also  be  used  to  analyze  the  drag  coefficient 
Cd-  Figure  3.10  shows  Cd  calculated  using  (3.2)  on  the  standard  2D  porous 
media  domain  for  a  variety  of  diameters.  This  result  resembles  a  Moody 
diagram  in  pipe  flow.  It  shows  linear  behavior  (as  previously  shown)  in  the 
low  Re  ranges  and  approaches  a  constant  at  high  Re  as  theory  suggests  should 
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(a)  Changes  in  Forchheimer  coefficient  for  Darcy-Forchheimer  equa¬ 
tion  from  least  squares  fit. 


Darcy-Forchheiemer  Fits 


(b)  Darcy-Forchheimer  plot  using  parameter  fit  over  different  ranges 
of  Re. 


Figure  3.7:  Darcy-Forchheimer  £t  results. 
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(a)  Changes  in  cubic  coefficient  for  cubic  equation  from  least  squares 
fit. 


(b)  Cubic  law  plot  using  parameter  fit  over  different  ranges  of  Re. 

Figure  3.8:  Cubic  law  fit  results. 
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Least  Squares  Fit  for  Power  Law  Coefficient 


Least  Squares  Fit  for  Power  Law  Power 


74 


(a)  Changes  in  power  term  coefhcient  for  (b)  Changes  in  power  for  power  equation 
power  equation  from  least  squares  fit.  from  least  squares  fit. 


Power  Law  Fits 


(c)  Power  law  plot  using  parameter  fit  over  different  ranges  of  Re. 


Figure  3.9:  Power  law  fit  results. 
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Figure  3.10:  Drag  coefficients  for  2D  porous  media  domains  with  circles  of 
varying  diameter. 

happen.  The  behavior  between  these  regions  is  interesting  and  requires  more 
study. 

It  is  important  to  show  that  our  numerical  method  converges  with 
spatial  rehnement.  Figure  3.11  shows  Cd  calculated  for  d  =  0.6  m  at  high  Re 
for  3  mesh  sizes.  This  plot  qualitatively  shows  that  our  calculated  values  of 
Cd  for  high  Re  do  converge  with  rehnement. 

3.5.5  3D  Emergent  Cylinder  Problem 

The  three-dimensionality  of  vegetation  has  a  major  effect  on  how  be¬ 
havior.  Flow  behaves  diherently  depending  on  whether  or  not  the  vegetation 
is  completely  submerged  or  whether  it  rises  above  the  huid  level.  If  it  is  com- 
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Figure  3.11:  Convergence  of  high  Reynolds  number  drag  coefficients  with  dif¬ 
ferent  mesh  sizes  h  for  2D  domain  with  d  =  0.6  m. 


Figure  3.12:  Surface  of  mesh  used  for  3D  calculations. 


56 


77 


(a)  Velocity  Magnitude 


(b)  Pressure  Contours 


(c)  u  Velocity  Contours 


(d)  Streamlines 


Figure  3.13:  Emergent  Cylinder  Re  =  0.025. 
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(a)  Velocity  Magnitude 


(b)  Pressure  Contours 


(c)  u  Velocity  Contours  (d)  Streamlines 


Figure  3.14:  Emergent  Cylinder  at  Re  =  60.15. 


58 


79 


Figure  3.15:  Darcy  and  non-Darcy  upscaling  results  for  3D  flow  around  a 
cylinder  of  diameter  0.7  m. 

pletely  submerged,  then  the  flow  height  is  important.  The  bed  also  has  a 
major  effect  on  the  flow. 

We  consider  a  cubic  domain  with  a  cylindrical  obstacle  as  seen  in  Figure 

3.2.  Figure  3.12  shows  an  example  of  a  tetrahedral  mesh  on  this  domain.  We 
incorporate  the  horizontal  periodic  boundary  conditions  as  described  in  Section 

3.3.  Figures  3.13  and  3.14  show  velocity,  pressure  contours,  u  velocity  contours, 
and  streamlines  for  flows  with  Re  =  0.025  and  Re  =  60.15  respectively.  In 
Figure  3.13  we  see  deformation  drag,  just  as  we  had  seen  in  2D.  In  Figure  3.14, 
we  qualitatively  see  surface  drag. 

We  also  consider  Darcy  and  non-Darcy  models  for  our  3D  results.  In  the 
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Figure  3.16:  Cd  over  a  range  of  Re  for  flow  around  3D  cylinders  of  different 
diameters. 

same  way  that  we  did  for  the  2D  results,  we  calculate  K  using  mathematical 
homogenization  and  perform  parameter  htting  to  estimate  parameters  for  non- 
Darcy  laws.  We  use  data  from  the  full  range  of  Reynolds  numbers  to  £t  for  the 
parameters  necessary  for  the  non-Darcy  laws.  The  resulting  drag  coefficients 
Cd  for  flow  around  a  cylinder  of  diameter  0.7  m  using  these  upscaling  models 
over  a  range  of  Re  are  shown  in  Figure  3.15.  Cd  over  a  large  range  of  Re 
for  different  diameter  cylinders  is  shown  in  Figure  3.16.  These  results  are 
qualitatively  similar  to  the  results  of  our  2D  analyses.  We  should  note  that 
at  these  drag  coefficients  are  substantially  different  in  magnitude,  but  not  in 
dynamics,  from  the  2D  results  at  low  i?e;  however,  those  domains  are  much 
less  porous  due  to  their  tight  packing,  causing  a  lower  theoretical  hydraulic 
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conductivity.  Darcy’s  Law  is  valid  for  very  small  Re  (<  1),  as  we  saw  in  2D. 
The  Darcy-Forchheimer  equation  and  the  power  law  (3.9)  match  the  Navier- 
Stokes  data  quite  well  over  the  whole  range  of  Re. 

These  3D  simulations  are  much  more  computationally  expensive  than 
their  2D  counterparts.  That  they  provide  similar  results  to  the  2D  simulations 
gives  credence  to  their  accuracy. 

3.5.6  3D  Submerged  Cone  Problem 

This  simulation  is  most  like  a  real  world  flow  over  submerged  vegetation. 
The  problem  uses  the  submerged  cone  domain  in  Figure  3.2.  It  models  the 
behavior  of  flow  through  a  periodic  bed  of  completely  submerged  plants.  The 
3D  behavior  of  this  flow  is  much  more  complex  than  that  in  the  previous 
example.  There  are  nontrivial  vertical  velocities  in  sections  of  the  domain.  By 
driving  flow  with  range  of  hydraulic  gradients  F  to  steady-state  we  see  the 
flow  structure  and  upscaled  results  for  a  variety  of  Reynolds  numbers.  The 
cones  in  these  simulations  have  a  height  of  0.8  m  and  base  radius  of  0.25  m 
and  are  located  in  the  center  of  the  cell.  For  our  scaling  parameter  d  we  use 
0.35  m,  the  diameter  at  the  center  of  mass. 

Figure  3.17  shows  results  from  a  low  Reynolds  number  flow  with  Re  = 
0.014.  The  length  scale  used  to  calculate  Re  is  the  diameter  of  the  base  of  the 
cone  which  is  0.4  m  in  this  case.  Notice  the  3D  behavior  in  the  streamlines. 
There  is  vast  flow  deformation  in  the  horizontal  and  vertical  directions  around 
and  over  the  cone.  In  the  vertical  slice,  we  see  a  thick  boundary  layer  around 
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the  bottom  and  the  cone  with  a  gradual  increase  in  velocity.  As  we  have  seen 
in  other  2D  and  3D  examples,  there  is  high  pressure  upstream  of  the  obstacle 
and  low  pressure  downstream.  Figure  3.18  shows  the  same  results  for  a  flow 
with  Re  =  6.3.  In  this  instance  the  deformation  is  less  prominent  farther 
away  from  the  cone  but  very  drastic  near  the  cone.  We  see  a  much  smaller 
boundary  layer  near  the  bottom  and  near  the  cone.  Above  this,  we  see  a 
standard  logarithmic  velocity  prohle.  This  displays  the  effects  of  form  drag  in 
the  vertical  direction. 

As  with  other  2D  and  3D  problems  done  here,  we  look  at  Darcy  and 
non-Darcy  upscaling  results.  In  Figure  3.19  we  see  similar  results  to  the  2D 
problems  and  the  3D  cylinder  problem.  At  low  Re  Darcy’s  Law  and  Darcy- 
Forchheimer  match  very  well  with  Navier-Stokes  data.  At  higher  Re,  Darcy’s 
Law  become  completely  invalid  and  Darcy-Forchheimer  and  a  power  law  pro¬ 
vide  reasonable  approximations  for  Cd-  The  full  range  of  simulations  were  used 
to  £t  for  the  various  non-Darcy  model  parameters. 

3.6  Conclusions 

Flow  through  porous  structures  and  vegetated  domains  can  be  quite 
complex.  For  very  low  Re  flows  Darcy’s  Law  accurately  and  effectively  models 
the  mean  flow.  Using  mathematical  homogenization,  the  resulting  hydraulic 
conductivities  can  be  relatively  easily  calculated.  However,  for  higher  Re  flow, 
it  is  not  as  simple.  Non-Darcy  models  such  as  Darcy-Forchheimer  provide  more 
accurate  results  than  Darcy’s  Law;  however,  they  require  important  parame- 
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(a)  Velocity  Magnitude 


(b)  Streamlines 


(c)  Pressure  Contours  (bottom  view)  (d)  Velocity  in  Vertical  Slice 

Figure  3.17:  Flow  characteristics  for  flow  around  a  submerged  cone  of  base 
diameter  of  0.5  m  and  a  height  of  0.8  m  at  Re  =  0.014. 
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(a)  Velocity  Magnitude 


(b)  Streamlines 


(c)  Pressure  Contours  (bottom  view)  (d)  Velocity  in  Vertical  Slice 

Figure  3.18:  Flow  characteristics  for  flow  around  a  submerged  cone  of  base 
diameter  of  0.5  m  and  a  height  of  0.8  m  at  Re  =  6.3. 
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Figure  3.19:  Darcy  and  non-Darcy  upscaling  results  for  3D  flow  around  a 
submerged  cone  with  base  diameter  0.5  m  and  a  height  of  0.8  m. 
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ters  to  be  estimated  and  are  not  effective  over  large  ranges  of  Re.  The  most 
descriptive  way  of  displaying  the  drag  effects  on  these  types  of  domains  over  a 
large  range  of  Re  is  a  chart  similar  to  a  Moody  diagram.  Creating  this  would 
require  performing  extensive  computational  simulations  over  many  different 
packings  of  vegetation  and  many  Reynolds  numbers.  This  approach  might  al¬ 
low  for  a  semi-automated  system  for  taking  remotely  sensed  geometry  informa¬ 
tion  (i.e.  using  LiDAR)  and  calculating  parametrized  resistance  factors.  Fur¬ 
ther  work  remains  for  issues  such  as  depth  dependance  with  submerged  vegeta¬ 
tion,  flexible  stems,  foliage,  and  potential  issues  like  wave/current /vegetation 
interaction.  Once  adequate  drag  analyses  are  performed,  the  upscaled  results 
could  be  used  to  model  flow  through  large,  complex  vegetated  domains  with¬ 
out  having  to  capture  the  hue  details  of  the  domain,  providing  acceptable 
estimations  of  flow  characteristics  while  being  computationally  viable. 
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Chapter  4 

Immersed  Structure  Approach 

4.1  Introduction 

The  work  presented  in  the  previous  chapter  describes  a  novel  method 
for  parameterizing  drag  for  flow  through  a  relatively  simple  domain  containing 
relatively  few  rigid  obstacles.  For  larger  scale  problems,  it  becomes  compu¬ 
tationally  unreasonable  to  resolve  all  of  the  complex  geometry  of  a  bed  of 
vegetation.  The  level  of  mesh  rehnement  necessary  to  accurately  resolve  on 
a  single  mesh  thousands  or  even  hundreds  of  pieces  of  vegetation  would  push 
the  limits  of  current  computational  capabilities.  Also,  as  mentioned  previously, 
plant  flexibility  is  often  an  important  factor  in  the  amount  of  resistance  to  flow 
due  to  vegetation.  To  incorporate  plant  flexibility  into  the  single-mesh  method 
introduced  in  the  previous  chapter  would  require  moving  meshes  and  3D  struc¬ 
ture  models  to  model  the  elastic  motion  of  the  vegetation.  This  would  require 
a  3D  single-mesh  fluid-structure  interaction  (FSI)  method  that  would  be  pro¬ 
hibitively  computationally  expensive  and  possibly  impossible  to  implement  at 
the  scale  that  we  require.  Instead,  in  this  chapter  we  present  an  immersed 
structure  approach  to  modeling  flow  over  larger-scale,  more  complex,  rigid 
and  flexible  beds  of  vegetation. 
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4.1.1  Bulk  Drag 


In  Section  3.2,  we  gave  a  general  overview  of  the  concept  of  drag  and 
important  concepts  like  the  drag  equation,  the  drag  coefficient,  and  how  it 
changes  as  Reynolds  number  changes.  The  drag  formulation  described  there 
is  for  ideal  flow  conditions.  For  a  larger-scale  vegetated  channel,  the  drag  de¬ 
pends  on  free-surface  effects,  turbulence,  and  complex  velocity  prohles.  These 
effects  are  described  in  detail  by  Petryk  [79].  Also,  the  effect  of  nearby  veg¬ 
etative  obstacles  affects  the  drag  pieces  of  vegetation.  As  described  by  Nepf 
[73],  for  a  vegetated  channel,  drag  force  per  unit  fluid  volume  is  dehned  by  the 
bulk  drag  equation 

Fdrag  =  ^pCdaU^  (4.1) 

where  Cd  is  a  non-dimensional  bulk  drag  coefficient,  p  is  the  fluid  density,  U 
is  the  mean  velocity,  and  a  is  the  projected  plant  area  per  unit  volume,  the 
so-called  vegetation  population  density.  Modeling  the  plants  as  cylinders,  a  is 
defined  as 

where  is  the  number  of  plants  per  unit  horizontal  area,  b  is  the  mean 
cylinder  diameter,  and  AS*  is  the  mean  spacing  between  plants,  and  H  is  the 
channel  depth. 

Whereas  for  a  single  obstacle  element,  the  drag  coefficient  Cd  is  only 
a  function  of  obstacle  shape  and  Reynolds  number,  the  bulk  drag  drag  coeffi¬ 
cient  is  also  a  function  of  population  density.  The  effect  of  population  density 
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can  even  be  observed  among  a  pair  of  cylinders  with  one  situated  downstream 
from  the  other.  Early  reviews  of  interference  between  two  cylinders  in  different 
arrangements  were  done  by  Zdravkovich  [116]  and  Zdravkovich  and  Pridden 
[118]  using  wind  tunnel  experiments.  It  has  been  observed  that  the  wake 
caused  by  the  upstream  cylinder  decreases  the  drag  coefficient  for  the  down¬ 
stream  cylinder.  This  has  been  shown  for  semi-inhnite  cylinders  by  Bokaian 
and  Geoola  [8]  and  Blevins  [7]  and  for  hnite-length  cylinders  by  Luo  et  al. 
[64].  This  decrease  in  drag  coefficient  is  caused  by  two  factors  of  the  wake. 
The  hrst  factor  is  that  the  downstream  cylinder  experiences  lower  fluid  veloc¬ 
ity  due  to  the  drag  of  the  upstream  cylinder.  The  second  factor  is  that  the 
wake  of  the  upstream  cylinder  results  in  lower  pressure  differential  around  the 
downstream  cylinder  according  to  Zukausras  [119]  and  Luo  et  al.  [64].  The 
decrease  in  downstream  drag  is  called  the  “sheltering  effect”  by  Raupach  [84]. 

There  is  a  long  history  of  attempts  at  calculating  the  bulk  drag  co¬ 
efficient  for  vegetated  channels.  Li  and  Shen  [56]  analyzed  the  experimental 
results  from  rigid  cylinders  due  to  Petryk  [79]  using  a  wake  superposition 
model.  They  found  that  further  downstream  in  the  array  of  cylinders  that 
the  bulk  drag  coefficient  asymptotically  approached  a  value  between  1.1  and 
1.2  at  a  Reynolds  number  of  9  x  10^  and  that  it  generally  slightly  decreased 
with  greater  spacing  between  the  cylinders.  Klaassen  and  Van  Der  Zwaard 
[48]  calculated  the  bulk  drag  coefficient  using  a  Chezy  formula  to  analyze  data 
for  air  flow  though  canopies  of  fruit  trees.  Their  estimates  for  Cd  are  higher 
than  1.2,  perhaps  due  to  more  turbulence  in  tree  canopies  than  cylinders  sub- 
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merged  in  water.  They  also  saw  the  bulk  drag  coefficient  decrease  with  greater 
spacing  between  obstacles.  Seginer  et  al.  [94]  performed  wind  tunnel  experi¬ 
ments  using  rigid  aluminum  cylinders  and  utilized  the  horizontal  momentum 
balance  equation  to  estimate  the  bulk  drag  coefficient.  They  concluded  that 
increased  obstacle  population  density  increased  the  turbulence  intensity  and 
thus  a  decrease  in  drag.  These  results  are  in  opposition  to  experimental  re¬ 
sults  by  Graf  and  Ko  [33]  and  the  earlier  discussed  results.  The  applicability 
of  these  hndings  have  been  questioned  by  Dunn  et  al.  [23].  Reid  and  Whitaker 
[88]  performed  experiments  with  rigid  cylindrical  obstructions  submerged  in 
channels  of  variable  depth.  They  used  Manning’s  equation  to  estimate  the 
bulk  drag  coefficient  for  different  depths,  estimating  Cd  to  be  1.77.  Burke  and 
Stolzenbach  [11],  using  a  RANS-type  computational  model  with  a.  k  —  e  tur¬ 
bulence  model  estimated  that  Cd  =  2.5  would  be  appropriate  for  all  velocities 
for  Spartina  grass  in  a  marsh.  They  mention  the  limitations  of  such  a  simple 
model  and  suggest  more  research  be  done  in  the  area.  From  experiments  using 
flexible  plastic  cylinders  Saowapon  and  Kouwen  [92]  estimated  the  bulk  drag 
coefficient  as 

C'd  =  2sin^0  (4.3) 

where  0  the  the  bent  angle  of  the  vegetation.  This  equation  which  is  given 
by  Hoerner  [36]  obviously  returns  values  between  0  and  2.  Den  Hartog  and 
Shaw  [22]  performed  experiments  over  a  flexible  corn  canopy  and  estimated 
mean  bulk  drag  coefficients  between  0.2  and  0.043.  They  provide  an  equation 
for  the  vertical  prohle  of  Cd-  These  early  studies  show  a  wide  variation  in 


70 


91 


calculated  values  of  Cd-  This  is  due  to  the  different  methods  of  calculating 
Cd,  experimental  setups,  and  measuring  devices.  While  they  do  not  provide 
concrete  values  or  even  ranges  of  Cd  for  specihc  flows,  the  important  consensus 
that  can  be  drawn  from  these  studies  is  that  Cd  generally  decreases  with  greater 
vegetation  spacing. 

4.1.2  Experimental  Methods  for  Calculating  the  Bulk  Drag  Coef¬ 
ficient 

For  verifying  our  methods,  we  focus  on  three  more  recent  very  well- 
respected  experimental  studies.  The  hrst  is  by  Dunn  et  ah  [23]  in  which  rigid 
and  flexible  cylinders  were  placed  in  a  flume  with  a  variety  of  spacings  and 
for  a  large  range  of  Reynolds  numbers.  They  used  acoustic  Doppler  velocime- 
ters  that  measure  velocity  and  turbulence  characteristics  which  they  used  to 
calculate  the  bulk  drag  coefficient.  The  second  study  is  by  Nepf  [73]  in  which 
wooden  cylinders  were  placed  in  a  variety  of  arrangements  in  a  flume  and 
acoustic  Doppler  and  laser  Doppler  velocimeters  were  used  to  measure  veloc¬ 
ities.  The  horizontal  momentum  balance  is  used  to  calculate  Cd-  The  third 
study  is  by  Wu  et  al.  [115]  and  uses  a  laboratory  wave  tank  to  measure  the 
wave  attenuation  by  vegetation.  They  perform  experiments  with  rigid  dowels, 
flexible  rubber  hoses,  and  actual  marsh  grasses. 

Using  Reynolds  averaging,  for  fully  developed  uniform  shear  flow  through 
cylindrical  obstructions,  Dunn  et  al.  [23]  take  the  horizontal  momentum  equa- 
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tion  to  be 

d  f  du  d  f  du 

°  j  +  "  j 

d  f  du  —r-;\  ,A  A\ 

+  &  “  j  • 

where  ”  denotes  the  time-averaged  value  and  '  denotes  the  fluctuating  compo¬ 
nent.  The  terms  pu'v' ^  and  pu'w'  are  called  the  Reynolds  stresses.  Using 

the  horizontal  averaging  technique  from  Raupach  and  Shaw  [86] ,  from  the  rate 


of  change  of  the  shear  stress,  we  get 


dz 


9S  -  \acyi 


(4.5) 


where  S  is  the  bed  slope,  Uh  and  Wh  are  the  velocities  averaged  in  the  horizontal 


plane,  and  is  the  horizontally  averaged  drag  coefficient.  Now,  solving  (4.4) 
for  U^,  we  have 


^  9S  -  ^.ju'w'f^) 
"  a/2uK^  ■ 


(4.6) 


They  calculate  two  versions  of  the  bulk  drag  coefficient.  One  uses  the  channel 


depth  H  as  the  vertical  length  scale  and  one  uses  the  mean  vegetation  height 


hv  as  the  vertical  length  scale.  These  bulk  drag  coefficients  are  calculated  by 

^  /o^  C'Xdz 

^  f^u^dz  ^ 

Jo 

and 


C,.  = 


Jo  ^  C'^u'^dz 

Jo"  uldz 


(4.8) 


Nepf  [73]  uses  a  different  horizontal  force  balance  relationship  to  calcu¬ 
late  Cd-  The  equation  is 

1  ~  W  f)M 

(1  -  ad)CBU^  +  -Cdad{-)U^  =  gh—  (4.9) 
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where  Cb  =  0.001  is  the  bed  drag  coefficient  and  ^  is  the  surface  slope.  The 
hrst  term  on  the  left  hand  side  is  the  is  bed  drag  and  the  second  term  is  the 
drag  contributed  by  the  stems. 

Wu  et  ah  [115]  calculate  the  bulk  drag  coefficient  by  analyzing  the 
wave  attenuation  due  to  vegetative  resistance.  The  theory  of  dissipation  of 
wave  energy  due  to  resistance  effects  has  been  developed  in  detail  by  Putnam 
and  Johnson  [81],  Bretschneider  and  Reid  [10],  and  Sleath  [96].  According  to 
Dalrymple  et  ah  [17],  assuming  linear  wave  theory,  the  conservation  of  wave 
energy  is 


d(Ec, 

dx 


=  -e,; 


(4.10) 


where  E  =  \pgA^  is  the  wave  energy  per  unit  area,  where  A  is  the  wave 
amplitude,  Cg  is  the  wave  group  velocity,  and  e„  is  the  dissipation  due  to  the 
vegetation.  Let  L  be  the  wavelength,  k  =  j L  be  the  wave  number,  and  H 
be  the  mean  water  depth.  Dalrymple  et  al.  [17]  and  Kobayashi  et  al.  [49] 
express  e„  in  terms  of  the  drag  force  from  the  vegetation: 


671  — 


N^{F^u  +  FyV  +  Fzw)dz, 


(4.11) 


Jo 

where  F^,  Fy,  F^  are  the  drag  forces  in  the  corresponding  directions  and  the 
overline  denotes  time-averaging.  Therefore,  we  have 

ev=  [  lrCdpN^by^\v\dz.  (4.12) 

Jo  ^ 

Mendez  and  Losada  [68]  assume  that  the  wave  heights  have  an  invariant 
Rayleigh  distribution  and  assuming  analytical  formulas  from  linear  wave  the¬ 
ory  obtained  the  formula  for  wave  energy  dissipation  per  unit  horizontal  area 


73 


94 


due  to  vegetation 


2uj 


sinh^  khy  +  3  sinh  kh^ 


^rp 


(4.13) 


?>k  cosh^  kh 

where  uj  is  the  wave  angular  frequency  and  Arms  is  the  root-mean-squared 
wave  amplitude.  For  a  single  wave  according  to  Dalrymple  et  al.  [17]  the 
solution  to  the  wave  energy  equation  (4.10)  gives  the  relationship 


A{x)  1 
Ai  1  +  ax 


(4.14) 


where  x  is  the  horizontal  coordinate  originating  at  the  beginning  of  the  vege¬ 
tated  area,  A^  is  the  incident  wave  height,  A{x)  is  the  wave  height  at  x,  and 
CK  is  a  damping  factor  dehned  by 

8  A  1  nr  1  sinh^  khy  +  Ssinhkhr 

=  — (^dAiONyk- — ; — — - — — - ,  ,  , 

Ovr  (smh  2kh  -|-  2kh)  smh  kh 

where  a  =  hr/H  for  submerged  vegetation  and  a  =  1  for  emergent  vegetation. 
Similarly,  due  to  Mendez  and  Losada  [68],  for  irregular  waves 


Ar 


X 


Ar 


1  -|-  da; 


(4.16) 


where 

2  ~  hAT  h  3 

^  {sinh  2kh  +  2kh)  sinh  kh 

Hence,  by  measuring  the  incident  wave  heights  and  wave  heights  at  various 

locations  within  a  vegetated  area,  the  bulk  drag  coefficient  may  be  calculated 

since  all  other  values  are  known  in  the  previous  equations. 
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4.1.3  Previous  Computational  Models  For  Flexible  Vegetation 

The  common  approach  for  computationally  modeling  flexible  vegetation 
is  by  treating  a  piece  of  vegetation  as  a  flexible  cantilever  beam.  Kutija  and 
Hong  [53]  first  proposed  this  method  using  standard  Timoshenko  [102]  beam 
theory.  Saowapon  and  Kouwen  [92]  developed  a  similar  model.  Darby  [18] 
describes  a  one-dimensional  model  that  may  be  used  for  both  flexible  and 
rigid  vegetation.  Erduran  and  Kutija  [25]  extend  the  work  of  Kutija  and  Hong 
[53]  by  using  a  3D  RANS  technique  with  a  combination  of  mixing  length  and 
eddy  viscosity  closure  schemes.  They  propose  quasi-3D  coupling  with  the 
shallow  water  equations.  Ikeda  et  al.  [40]  use  cantilever  beam  theory  and 
2D  LES  to  model  “monami,”  the  waving  of  flexible  plants  due  to  large  eddies 
caused  by  instabilities  in  the  flow  held.  They  utilize  a  separate  vegetation 
grid  to  track  the  movement  of  each  piece  of  vegetation.  Monami  is  explained 
and  modeled  in  depth  by  Nepf  and  Ghialberti  [72].  Velasco  et  al.  [106]  use 
classical  beam  theory  to  compute  the  displacement  of  hexible  beams  under 
small  to  moderate  dehection.  Kubrak  et  al.  [52]  extend  the  method  developed 
by  Kutija  and  Hong  [53]  for  large  dehections.  Li  an  Xie  [54]  build  off  work 
involving  stiff  vegetation  by  Li  and  Zeng  [55],  add  a  thin  plate  of  “foliage”  to 
the  stems  and  use  a  3D  LES  scheme  with  a  Smagorinsky  closure  model.  They 
do  not  completely  couple  these  models.  Instead,  they  used  a  hnite  difference 
method  developed  by  Al-Sadder  and  Al-Rawi  [1]  to  solve  the  beam  equation 
for  a  large  variety  of  inhow  velocities  and  performed  a  parameterization  for 
the  bent  stem  height  based  on  how  velocity. 
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All  of  the  above  methods  only  consider  planar  bending;  however  in  a 
flow  with  large  velocity  or  a  great  deal  of  directional  change  or  turbulence, 
we  would  not  expect  planar  bending.  Also,  most  of  these  models  do  not 
appropriately  handle  large  deflections  and  are  only  valid  for  small  deflections. 
They  also  tend  to  not  be  tightly  coupled  with  the  flow  solver.  We  want  to 
develop  a  method  using  cantilever  beam  theory  to  model  flexible  vegetation 
that  is  as  robust  and  accurate  as  possible  but  is  not  overly  computationally 
expensive. 

4.2  Elastic  Cantilever  Beam  Theory 

4.2.1  Description  of  Problem 

We  model  flexible  vegetation  as  inextensible  cantilever  beams,  beams 
that  are  anchored  on  exactly  one  end.  Most  theory  for  elastic  cantilever  beams 
is  only  valid  for  small  deflections;  however,  the  deflection  of  vegetation  due  to 
the  force  of  a  fluid  flow  is  usually  relatively  large.  Therefore,  we  must  rely  on 
more  general  beam  theory  which  is  valid  for  large  deflections.  The  classical 
theory  of  the  elastica  is  described  by  Love  [62] .  Early  analytical  and  numerical 
work  involving  the  large  deflection  of  cantilever  beams  was  done  by  Bisshopp 
[6].  Deriving  from  this  classical  theory  Steigmann  and  Faulkner  [97]  outline 
the  variational  theory  for  inextensible  cantilever  beams.  This  is  general  theory 
that  is  valid  for  small  and  large  deflections. 

We  describe  the  3D  cantilever  beam  in  terms  of  an  arclength  parameter 
s  G  [0,  L],  where  L  is  the  length  of  the  beam.  We  dehne  an  orthonormal  basis 
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Figure  4.1:  An  inextensible  cantilever  beam. 


{ei(s),  62(5),  63(5)}  which  describes  the  orientation  of  the  beam’s  center-line  at 
s  and  a  vector  r(s)  which  gives  the  position  of  the  center-line  at  s  in  physical 
space.  ei(s)  is  the  unit  vector  in  the  direction  tangent  to  the  center-line 
increasing  in  arclength.  62(5)  and  63(5)  span  the  cross  section  of  the  cylinder 
at  s.  We  refer  to  {ej(s)}  as  the  material  basis  and  to  r(s)  as  the  material 
location.  We  dehne  a  local  basis 

{e°}  =  {ei}s=o  (4.18) 

that  will  serve  as  a  hxed  basis  along  the  length  of  rod.  In  the  undeformed 
state,  r(s)  =  X(s)  and  {ej(s)}  =  {Ei(s)}.  Let  '  denote  differentiation  by  s. 
Since  the  beam  is  inextensible,  we  have 

r'(s)  =  ei(s).  Vs  e  [0,  L].  (4.19) 
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The  given  constitutive  conditions  assure  that  deformations  from  the  un¬ 
deformed  state  to  the  material  state  are  inextensible  and  orientation-preserving. 
We  assume  that  the  cross  section  of  the  beam  at  s  remains  planar,  suffers  no 
strain,  and  is  orthogonal  to  ei(s).  Since  {ej(s)}  and  {Ej(s)}  are  right-handed 
orthonormal  bases,  then  we  can  dehne  a  rotation  transformation 


R(s)  =  ej{s)  0  Ej(s), 


(4.20) 


such  that 

ej(s)  =  R(s)Ei(s).  (4.21) 

Let  Ki(s)  be  the  twist  per  unit  length,  and  let  ^2(5)  and  ^3(5)  be  the 
curvature  about  the  ei(s)  and  62(5)  axes  respectively.  Let  k{s)  =  fi;i(s)ei(s) -t- 
^2(5)62(5)  -h  ^3(5)63(5).  Then 


e[  =  K  X  Bj. 


Suppose  the  undeformed  conhguration  has  twist  and  curvatures  ^^(s),  ^^(s), 
and  ^3(5).  Note  that  these  are  not  necessarily  zero.  Similarly,  we  have  the 
relationship 


(4.22) 


where  e  is  the  permutation  symbol. 


Let  W (k,  kP,  s)  be  the  strain  energy  per  unit  arclength  as  a  function  of 
curvatures.  Similarly  let  U (e*,  e'-,  s)  be  the  strain  energy  per  unit  arclength  as  a 
function  of  e*  and  e'.  It  is  obvious  that  we  must  have  W {k,  s)  =  U{ei,  e'-,  s). 
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The  strain  energy  of  the  system  is  then 

i*L  pL 

S{ei)  =  Uds=  Wds.  (4.23) 

Jo  Jo 

Suppose  that  f  is  an  end  load  at  s  =  L  and  q  is  the  piecewise  continuous  dis¬ 
tributed  load  along  the  arclength  of  the  beam.  The  potential  energy  functional 
of  the  system  is 

E{r,ei)  =  S{ei)  -  f  q  ■  rds  -  f  ■  r(L).  (4.24) 

Jo 

States  {r,  e*}  that  minimize  E  are  static  states  for  the  cantilever  beam  under 
the  given  loads. 

It  can  be  shown  using  the  calculus  of  variations  that  the  problem  of 
hnding  stable  minima  for  E  is  equivalent  to  solving  the  Euler-Lagrange  equa¬ 
tions 

F'(s)  +  q(s)  =  0  (4.25) 

M'(s)  =  F(s)  X  ei(s)  (4.26) 

with  the  boundary  conditions  F(L)  =  f  and  M(L)  =  0,  where  M  = 
and  F  is  a  vector  of  Lagrange  multipliers.  We  notice  that  these  equations 
are  exactly  the  equations  of  equilibrium  of  forces  and  equilibrium  of  moments 
respectively,  where  F  is  the  force  vector  and  M  is  the  vector  of  moments. 

The  most  important  theory  for  the  bending  of  beams  is  the  Euler- 
Bernoulli  Law,  which  linearly  relates  the  curvature  and  twist  of  the  beam  to 
the  bending  and  twisting  moments  respectively.  The  Euler-Bernoulli  Law  is 

M  =  GJ{^Ki  —  K^)®!  T  EIi[k2  —  ^2)^2  +  —  ^3)63,  (4.27) 
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where  GJ  is  called  the  torsional  rigidity  and  EI2  and  EI^  are  called  the  flexural 
rigidities  about  62  and  e^.  E  is  Young’s  modulus  and  I2  and  I2  are  the  area 
moments  of  inertia  about  62  and  ei  respectively.  G  is  the  shear  modulus  and  J 
is  the  polar  second  moment  of  inertia  about  ei.  Assuming  the  Euler-Bernoulli 
Law,  we  now  have  the  strain  energy  per  arclength  as 

W=^  [GJ{ki  -  +  Eh{K2  -  4?  +  .  (4.28) 

We  desire  to  find  the  material  basis  and  location  under  a  given  dis¬ 
tributed  load.  For  our  purposes,  assume  that  q  is  a  known  function  of  x,  y 
and  and  that  f  =  0.  We  describe  and  analyze  a  variety  of  numerical  methods 
for  this  problem.  Some  methods  involve  numerically  solving  the  equilibrium 
equations  (4.25)  and  (4.26)  using  the  Euler-Bernoulli  Law  (4.27)  as  a  constitu¬ 
tive  equation.  Others  involve  finding  a  solution  that  minimizes  the  potential 
energy  functional  (4.24). 

4.3  Numerical  Techniques 

We  propose  a  variety  of  acceptable  techniques  for  solving  the  system 
described  in  the  previous  section.  All  of  the  previous  work  discussed  in  Section 
4.1  uses  models  that  are  not  nearly  as  general  as  this.  Hence,  they  tend  to  be 
unable  to  handle  large  deflections.  They  also  only  allow  planar  bending.  We 
want  to  allow  full  3D  bending.  Allowing  very  large  deflections  and  3D  bending 
will  give  us  as  robust  a  method  as  is  possible. 
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4.3.1  Nonlinear  Finite  Element  Method 


In  this  method,  we  use  a  nonlinear  finite  element  method  to  find  local 
extrema  of  the  potential  energy  functional  (4.24).  A  version  of  this  technique 
for  planar  bending  is  presented  by  Fried  [27].  Assume  that  the  distributed 
load  q  =  [o'!,  5'2,  Q's]^  is  a  known  function  of  spatial  coordinates  x,  y,  and  i:. 
Assume  that  there  are  no  end  loads.  We  define  three  Euler  angles  6,  0,  and  ip 
such  that 


L 

^2 


L 

L 


sin  6*  cos  0, 
sin  9  sin  ip, 
cos  9 


(4.29) 

(4.30) 

(4.31) 


where  9  =  0  corresponds  to  when  the  center-line  is  along  the  z-axis. 

We  can  write  the  twist  and  curvatures  in  terms  of  9,  (p,  and  ip  as 


/tl 

=  (p'  +  Ip'  COS  9, 

(4.32) 

«2 

=  9'  sin  (p  —  Ip'  sin  9  cos  0,  and 

(4.33) 

«3 

=  9'  cos  (p  +  Ip'  sin  9  sin  0. 

(4.34) 

Let  us  assume  that  the  undeformed  state  has  zero  curvatures  and  twist  for  all 
s.  The  potential  energy  of  the  system  (4.24)  becomes 

E{9,(p,ip)  =  [  [I- {G  J{(p' +  Ip' cos  9Y  +  EIi{9' sin  (p  —  Ip' sin  9  cos  (py 

Jo  2 

-|-  El2{9'  cos  (p  +  Ip'  sin  9  sin  0)^)  —  qiTi  —  g2i’2  —  (4.35) 
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Suppose  qi  =  Q[,  q2  =  Q2,  qs  =  Q's,  QiW  =  ^2(1)  =  Q3(l)  =  0,  and 
ri(0)  =  r2(0)  =  r3(0)  =  0.  Then  (4.35)  becomes 

=  [  [^{G  J{4>' +  Ip' cos  6)'^  +  EIi{6' sin  (p  —  Ip' sin  6  cos  4>)'^ 

Jo  2 

+  El2{6'  cos  (p  +  'ip'  sin  6  sin  0)^)  +  Qi  sin  6  cos  tp 
+  (^2  sin6*sin0  +  Qs  cos6*](is.  (4.36) 

For  simplicity,  assume  that  Eli  =  EI2  =  El,  which  would  be  the  case  for  a 
homogeneous,  isotropic  material  with  a  circular  cross  section.  This  gives 

/■i  1 

E{9,(p,'ip)  =  /  [-(GJ(0'^(1  +  cos^  6*)  +  2-0'0' cos6')  +  i?/(6''^  +  0'^  sin^  6*) 

Jo  2 

+  Qi  sin  6*  cos^  +  (^2  sin  6*sin0  +  Qa  cos6*](is.  (4.37) 

Using  quadratic  hnite  elements  we  calculate  the  element  potential 
energy  E^..  The  element  has  3  nodes  with  two  on  the  element  boundary  and 
one  at  the  element’s  midpoint.  The  reference  element  is  the  range  [—1, 1],  and 
the  shape  functions  on  the  reference  element  are 

w‘K)  =  |i«K-l).l-y.i?«  +  l)].  (4.38) 

Interpolation  of  the  Euler  angle  functions  on  the  element  is  expressed  by 

9  =  Oe^w,  tp  =  'ipe^w,  and  (p  =  (pj^w,  (4.39) 

where  0e^  =  (6^1,  6^2,  6^3)  where  6*1  =  6*(— 1),  6^2  =  6^(0),  and  9^  =  9{1),  is  the 
element  nodal  values  vector  and  likewise  for  tp  and  (p.  Let 

=  [Ol,  01,  01,  O2,  02,  02,  O3,  03,  03]^  (4-40) 
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be  the  element  Euler  angle  vector. 

Let  the  element  size  be  h,  so  ds  =  \hd^,  (  )'  =  2h~^{  ),  and  (  )  = 

We  integrate  using  2  point  Gauss  quadrature  with  Gauss  points,  Ci  =  —  1/v^, 
^2  =  1/v^-  Let  the  index  j  =  1,  2  be  such  that  9j  is  the  value  of  9  at  Cj,  so 

=  ^rw(Cj),  9j  =  0^w{Q.  (4.41) 


For  brevity  dehne  wj  =  "w{(j)  and  wj  =  w{(j).  The  approximate  element 
potential  energy  functional  is  now 
1  .  ^ , 

Ee  =  -  sin^  9j)  +  +  cos^  9j)  +  2ipj^j  cos  9j)].  (4.42) 

^  i=i 

We  dehne  the  element  gradient 


ge  = 


dE, 

d^e 


(4.43) 


and  the  element  stiffness  matrix 


dge  ^  d‘^Ee 
d^e  d^l  ' 


(4.44) 


We  want  to  solve  the  problem  g  =  0.  Our  global  initial  angle  vector  is 
$0-  Use  (4.43)  to  assemble  the  initial  gradient  vector  go.  Delete  the  entries 
that  correspond  to  hxed  points.  Assemble  the  initial  stiffness  matrix  Kq.  Use 
the  Newton-Raphson  method 


^Ti+l  Kn  gfi 


(4.45) 


until  the  residual  is  within  a  given  tolerance.  Since  this  is  an  optimization 
problem,  we  must  verify  that  we  are  indeed  at  a  minimum.  If  all  of  the 
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eigenvalues  for  Kjv  where  N  is  the  terminal  Newton  iteration  are  positive, 
then  we  are  at  a  local  minimum.  We  can  also  verify  that  (4.25)  and  (4.26)  are 
satished. 

Since  this  problem  is  highly  nonlinear,  it  is  possible  for  the  Newton- 
Raphson  method  to  diverge  or  never  converge  if  the  initial  guess  is  greatly 
different  than  the  solution.  Thus,  it  is  best  to  use  incremental  loading  if  the 
load  is  going  to  cause  a  large  deflection.  Results  using  incremental  loading  for 
a  variety  of  uniform  distributed  loads  are  shown  in  Figure  4.2.  Notice  the  3D 
behavior  of  the  bending.  The  loads  are  greater  in  the  y  direction  so  we  see  a 
greater  deflection.  This  is  the  type  of  behavior  that  cannot  be  captured  with 
planar  bending  methods. 

4.3.2  Linearized  Method 

Rather  than  solving  the  above  minimization  problem,  we  have  shown 
above  that  it  is  equivalent  to  solving  (4.25)  and  (4.26),  the  equilibrium  of  forces 
and  moments  respectively  with  boundary  values  M(0)  =  0  and  F(L)  =  F. 
(4.25)  can  be  solved  for  F  by  simple  numerical  integration  techniques  if  q  is 
known. 

Assume  that  the  beam’s  unloaded  state  is  straight,  i.e.  /t°(s)  =  0  V  a;  G 
[0,L]  and  that  =  [1,0,0]^,  e®  =  [0,1, 0]"^,  and  elj  =  [0,0,1]^.  Let  the 
curvatures  be  dehned  as  in  (4.32,4.33,4.34).  Let  ^*(s)  =  [-^(s),  6*(s),  0(5)]"^ 
and  ^*^(s)  =  [ip' {s) ,  9' {s) ,  (j)' {s)]'^ .  The  dehnitions  above  imply  that  $(0)  =  0. 
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Bending  in  x-z  plane  due  to  uniform  distributed  loads 


X 


Bending  in  y-z  plane  due  to  uniform  distributed  loads 


y 


Figure  4.2:  Deflections  in  the  x-z  and  y-z  planes  due  to  uniform  distributed 
loads.  This  has  been  solved  using  the  nonlinear  hnite  element  energy  mini¬ 
mization  method. 
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Using  (4.27),  (4.26)  becomes 

G J k\Q^-\-E I\K2G2~^^ ^ ~  R(Fxei)  (4.46) 
where 

cos  9  cos  (j)  cos  6  sin  0  —  sin  6 

R  =  sin  Ip  sin  6  cos  0  —  cos  ip  sin  cp  sin  ip  sin  6  sin  cp  +  cos  ip  cos  (p  cos  6  sin  ip 

cos  Ip  sin  9  cos  (p  +  sin  ip  sin  (p  cos  ip  sin  9  sin  (p  —  sin  ip  cos  (p  cos  9  cos  ip 

(4.47) 

is  the  rotation  matrix.  R  is  orthonormal,  so  R“^  =  R^.  Let 

■  GJ  0  o' 

E  =  0  EU  0  (4.48) 

0  0  EIz 

be  a  diagonal  matrix  with  the  stiffness  values  on  the  diagonal.  Let 

[H^,H2,H^Y  =  n{Yxe,). 

(4.46)  becomes 

{Ell  —  EIz)  1^3 

EiK  =  {EIz  ~  GJ)  Ki  Ks  +  Hs  .  (4.49) 

{GJ  —  Ell)  ^2  —  Ez 
Let  the  right-hand  side  be  known  as  H.  Let 

1  0  sin  6* 

G  =  0  cosip  cos9smip  .  (4.50) 

0  —  siii'^  cos  9  cos  Ip 

det(G)  =  cos 6*,  so  G  is  invertible  if  6*  7^  n|  for  n  =  ±1,  ±3, ...  .  Now  (4.49) 
becomes 

^G$'  =  E^^H.  (4.51) 

as 
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Using  the  product  rule,  we  have 

G$"  +  G'$'  =  E-^H  (4.52) 

Assuming  det(G)  7^  0,  this  can  be  written  as 

=  -G-^G'$'  +  G'^E-^H.  (4.53) 

The  condition  M(L)  =  0  implies  ^'{L)  =  0,  so  this  is  a  system  of  second 
order  initial  value  problems  that  is  solvable  with  standard  techniques  such 
as  Runge-Kutta  methods.  This  is  discussed  by  Raboud  et  ah  [82].  This 
method  is  fast  and  accurate  for  small  deformations;  however,  we  are  unable  to 
perform  incremental  loading  for  the  system  of  ODEs,  since  it  is  not  an  iterative 
procedure.  Since  the  system  is  so  nonlinear,  for  good  accuracy  for  problems 
with  large  loads,  it  requires  a  great  deal  of  segmental  steps  for  an  ODE  solver. 
By  adding  another  boundary  condition  $(0)  =  0,  which  implies  that  the  beam 
is  straight  at  the  base,  we  get  a  boundary  value  problem. 

When  there  are  large  deflections  and  we  need  to  use  incremental  loading 
to  get  an  accurate  result,  we  can  use  the  hnite  element  method.  Dehne  B  = 
(0,L).  Let  (•>  ')b  be  the  L2  inner  product  on  B.  Let 

S  =  {v|nj  G  C^{B)  for  i  =  1,2,3  and  v(0)  =  0  and  v'(L)  =  0}. 

Let  W  =  s.  We  call  $  G  S'  a  weak  solution  to  (4.53)  if 

($",  w)g  =  (^-G-^G'$'  +  G-^E-^H,  w)  ^  (4.54) 
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Vw  G  W.  Integrating  by  parts,  we  have 

-(#',wOs=  (^-G-iG'#'  +  G-^E-iH,w)^  (4.55) 

We  now  linearize  this  problem  so  that  Picard  iteration  may  be  used.  Lag  the 
values  in  the  right  hand  side  to  solve  for  the  new  $. 

-  =  (-G„-‘G:<I>;.  +  G;'E-‘H„,w)^  (4.56) 

If  the  sequence  {^*n}  ^  as  n  — >  cx)  then  $  is  a  weak  solution  to  the  original 
problem.  Suppose  that  EI2  =  EI3  =  EL 

Lemma  4.3.1.  //u  G  C^{B),  then 

<  L\u\hI{B)-  (4-57) 


Proof.  u(a;)  =  u'{s)ds,  so, 


px  pL 

l^('5)l  =  I  /  u'(s)(is|  <  /  |u'(s)|(is. 

Jo  Jo 


By  Cauchy-Schwarz, 


\u'{s)\ds  <  ||l|U2(s)||||u'||i2(B)  = 


Squaring  both  sides,  integrating  over  B,  and  taking  a  square  root  gives. 


||u||l2(b)  <  L|u|hi(b). 


□ 


£^2  I  JP  I 

Theorem  4.3.2.  If  is  sufficiently  small,  the  sequence  {$n}  converges. 
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Proof.  Define  the  smooth  nonlinear  map  h  ; 

h($)  =  +  G-^($)E-^H($).  (4.58) 

Now,  (4.56)  becomes 

-  =  (h($„),w)^.  (4.59) 

By  linearity, 

-  ($n+l  -  =  (h($„)  -h($„_l),w)g. 

Since  ^n+i  —  G  S'  =  W,  let  w  =  $n+i  —  to  get 

(sy,  -  «y,  -  «;)  ^  =  -  (h(4>„)  -  h(#„_i),  #„+,  -  . 

Now, 

=1$  , ,  —  $  P  ,  >  J_  11$  ,  1  —  $  ll^o 

t^n+l  ^nt^n+l  ^njj^  l^n+1  ^n\}ji  —  ^2  ^nWip.) 

by  Lemma  4.3.1.  So  now,  we  have 

^  ||$„+1  -  $n||5^2  <  |(h($n)  -  h($„_i),  ^n+1  “  ^n)Bl  • 

By  the  Cauchy-Schwarz  inequality, 

|(h($„)  -  h($„_i),  -  $„)g|  <  ||h($„)  -  h($„_i)||^2  ||$n+l  -  ^n\\L2  , 

so  now  we  have 

||$„+i  -  #,11^2  <  ||h($„)  -  h($,,_l)||^2  . 
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By  Taylor’s  Theorem,  there  exists  G  S  such  that 


h(#„)  =  h(<I>„_0  +  (  ^h(<l.*)  I  («„  - 


Therefore, 

It  can  be  shown  that 


d 


h($) 


|#n  -  ^n-l| 


L2 


L°° 


El  cos^  9 


<  C  =  0(1). 


Finally,  we  have 


L^IFIO 


l^n+l  0 


-  $n_l||^2  ,  (4.60) 

SO,  if  CO <  1  then  we  have  a  contraction.  By  the  Banach  hxed-point 
theorem,  the  sequence,  {^*n}  converges. 


□ 


Thus,  we  are  ensured  convergence  using  Picard  iteration  if  the  ratio  of 
the  beam  force  times  the  square  of  the  beam  length  and  the  beam  flexibility 
is  small  and  the  system  is  not  near  a  state  where  cos  6  =  0.  Notice  that 
these  angles  are  also  the  angles  that  cause  G  to  be  singular.  This  method  is 
especially  unstable  when  6  is  near  such  an  angle. 

Let  the  Galerkin  hnite  element  spaces  G  S  and  C  hF  be  dehned 
using  quadratic  hnite  elements  in  the  standard  way.  The  element  has  3  nodes 
with  two  on  the  element  boundary  and  one  at  the  element’s  midpoint.  The 
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reference  element  is  the  range  [—1, 1],  and  the  shape  functions  on  the  reference 
element  are 

vf  K)  =  K.  vl,  -  1),  1  -  +  1)],  (4.61) 

Let  $1,  $2  and  $3  designate  the  values  of  $  at  nodes  1,  2  and  3  respectively 
and  let  #e  =  The  Galerkiu 

problem  is  to  hnd  G  such  that 

-  +  G-^E-^H,  (4.62) 

Vw^  e  W^. 

Let  n  be  the  iteration  number.  Let  Ue  be  the  number  of  beam  elements 
and  Uq  be  the  number  of  quadrature  points  used  on  each  element.  Let  Cg  = 
[ci,C2,C3].  Let  So,  si, ...,  Sng  be  the  nodes  on  the  element  boundaries.  Let 
he  be  the  element  diameters.  Now,  suppose  Wg  =  CeVg.  Let  Uj  and  Q  for 
j  =  1,  ..,nq  be  the  quadrature  weights  and  points  for  an  Uq  point  quadrature 
rule.  Now  we  have 


0  =  (-G-'G'#li  +  G-iE-iH,w^j^ 

rie  —  l 

=  E 

i=0 


$e,nVe  '  CgV^ 


-G-lG'$e,n-lVe^'  +  G-'E-'h)  ■  CeV^jds 


rie  —  l 


i=0  j  =  l 

h 


■  CeV* 


-G~^G'#e,n-lV^'  +  G-^E-^H)  ■  CeV^]|,=,;.  (4.63) 
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This  is  only  true  for  all  Ce  if  for  each  element, 


0  = 


+ 


(4.64) 


Let 


de  =  [di,d2,d3]  =  ^ 


i=i 


h, 


\s=Cj 


(4.65) 


and 


Ke  —  7 

i=i 


—  v'^Vv^')'^ 


L=Cj 


(4.66) 


Ke  is  the  element  stiffness  matrix,  and  de  is  the  matrix  of  element  load  vectors. 
Let  Oe  =  [6'i,6'2,6'3]'^,  and  =  [01,02,03]^- 

Now,  form  the  global  stiffness  matrix  K,  load  vectors  di,  d2,  and  d3 
and  solution  vectors  V’,  and  (p.  Impose  the  boundary  conditions  by  deleting 
the  first  row  and  column  of  K  and  hrst  entry  of  the  load  and  solution  vectors. 
This  gives  us  the  matrix  problems 


K'lP  =  di. 

(4.67) 

K0  =  d2,  and 

(4.68) 

KcP  =  ds. 

(4.69) 

Notice  that  these  3  problems  are  completely  decoupled.  In  the  full  nonlinear 
hnite  element  problem,  the  linear  system  we  had  to  solve  was  3n  x  3n  where 
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n  is  the  number  of  nodes.  Here  we  solve  3  n  x  n  systems,  and  we  never 
have  to  update  the  stiffness  matrix.  While  the  linearized  method  may  take 
more  Picard  iterations  than  the  nonlinear  hnite  element  method  takes  Newton 
iterations,  each  iteration  for  the  linearized  method  is  much  less  expensive. 


Using  a  2  point  Gauss  quadrature  rule,  with  cui  =  a;2  =  1,  Ci 
and  C2  =  1/ Vs, 

7  -8  1  ' 

-8  16  -8  , 

1  -8  7 


-1/V3 


(4.70) 


4.3.3  Comparison  of  Methods 

As  shown  in  4.2.1,  the  systems  solved  by  the  fully  nonlinear  method 
and  the  linearized  methods  presented  above  are  equivalent.  However,  there  are 
major  differences  with  the  numerical  methods  presented.  For  each  iteration, 
the  fully  nonlinear  method  is  much  more  expensive.  The  stiffness  matrix  K 
must  be  assembled  for  each  iteration  and  a  3n  x  3n,  where  n  is  the  number  of 
nodes,  linear  system  must  be  solved.  For  the  linearized  method,  the  problem 
can  be  decoupled  into  3  independent  problems  for  each  step.  Also,  the  stiffness 
matrices  are  the  same  for  each  step,  so  they  only  need  to  be  formed  once. 
This  is  also  benehcial  for  the  linear  solve,  as  an  LU  decomposition  may  be 
performed  at  a  cost  of  ~  floating  point  operations  which  allows  the  system 
to  be  solved  by  back  substitution  at  a  cost  of  ~  floating  point  operations. 
Therefore  the  total  floating  point  operations  for  a  Picard  iteration  for  the 
linearized  method  is  ~  6n^,  where  it  is  ~  18n^  for  a  Newton-Raphson  iteration 
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for  the  nonlinear  method.  This  excludes  the  cost  of  forming  the  stiffness  matrix 
for  the  nonlinear  method,  which  could  be  quite  high.  The  approximations  for 
counts  of  floating  point  operations  for  linear  algebra  are  due  to  Trefethen  and 
Ban  [103] .  It  should  be  noted  that  for  the  nonlinear  and  linearized  systems  that 
the  stiffness  matrices  are  sparse  banded  matrices,  for  which  there  are  many 
iterative  methods  one  would  use.  However,  we  will  generally  be  using  relatively 
few  beam  elements,  so  our  stiffness  matrices  will  be  of  small  dimension  where 
dense  direct  solver  techniques  are  less  expensive  than  iterative  ones.  If  one 
wanted  to  model  beams  with  a  large  number  of  elements,  iterative  solvers 
would  surely  be  a  good  choice;  however,  we  generally  see  good  results  with  a 
small  number  of  beam  elements. 

Both  methods  have  a  signihcant  problem  being  of  unstable  if  0  if 
n  =  ±1,  3.  When  9  has  that  value  somewhere  along  the  beam  for  the  nonlinear 
method,  the  stiffness  matrix  K  becomes  singular.  For  the  linearized  method, 
it  forces  division  by  zero  on  the  right- hand- side.  Therefore,  we  must  check  to 
see  if  9  is  close  to  one  of  these  values  before  trying  to  solve  the  system.  If 
it  is  close,  then  we  add  a  small  amount  of  noise  to  9  to  avoid  the  instability. 
The  Newton-Raphson  and  Picard  iterations  methods  tend  to  be  forgiving,  and 
usually  adding  the  noise  does  not  hinder  convergence. 

A  useful  property  of  both  methods  is  that  they  are  designed  to  handle 
incremental  loading.  For  large  loads,  if  the  load  is  applied  to  the  system  and  the 
initial  guess  is  not  close  to  the  stable  equilibrium,  the  methods  are  very  likely 
to  diverge.  However,  if  the  load  is  incrementally  increased  and  the  equilibrium 
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state  from  the  previous  increment  is  used  as  an  initial  guess,  then  it  is  highly 
likely  that  the  nonlinear  iteration  will  converge.  Taking  incremental  steps  also 
may  cause  the  method  to  converge  with  fewer  total  nonlinear  steps.  Each 
increment  is  likely  to  converge  in  very  few  steps;  whereas,  with  no  incremental 
loading  the  system  may  jump  around  for  many  steps  before  converging. 

4.3.4  Beam  Benchmark  Problems 

We  test  our  two  solution  methods  for  large  deflections  of  elastic  can¬ 
tilever  beams  with  several  benchmark  problems.  These  problems  test  the  abil¬ 
ity  of  our  methods  to  handle  very  large  deflections  due  to  end  loads,  uniform 
and  variable  distributed  loads,  and  mixtures  of  types  of  loads.  It  is  important 
that  we  verify  that  they  provide  accurate  results  for  three-dimensional  bend¬ 
ing.  We  do  this  by  testing  all  bending  in  the  x-z  plane,  the  y-z  plane,  and  the 
plane  45  degrees  between  the  two.  We  are  also  interested  in  convergence  rates 
and  the  number  of  nonlinear  iterations  required  for  convergence. 

It  is  standard  for  the  test  problems  to  non-dimensionalize  the  equations. 
We  do  that,  in  effect,  by  taking  the  length  of  the  beam  L  to  be  1  m  and  the 
flexibility  El  to  be  1  Nm^  in  our  numerical  methods.  The  test  problems  give 
non-dimensional  end  loads  F  and  distributed  loads  q.  We  take  the  nonlinear 
tolerance  used  to  test  convergence  to  be  1.0  x  10“^  for  all  cases.  The  initial 
guess  for  all  methods  will  be  an  undeformed,  straight  beam. 

The  hrst  benchmark  problems  that  we  look  at  are  from  Fried  [27] .  Their 
method  is  well-validated  for  large  deflections  and  is  similar  to  a  2D  version  of 
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the  nonlinear  method  presented  here.  The  first  test  takes  Fi  =  1.5  and  all  other 
loads  to  be  0.  They  use  4  cubic  elements,  and  we  use  4  quadratic  elements. 
Table  4.1  shows  the  end  deflections  for  the  converged  system  for  the  hrst  6 
iterations  for  Fried’s  [27]  method  and  our  two  methods.  We  see  that  after  only 
a  few  iterations  both  of  our  methods  are  extremely  close  to  their  converged 
values  and  they  match  up  very  well  with  the  previous  results.  The  converged 
deflections  are  0.4109928  m  for  the  previous  results  and  0.41097726  m  for  both 
of  our  methods.  It  is  important  to  note  that  this  problem  is  a  large-deflection 
problem  and  one  that  small-deflection  methods  would  generally  fail  at  solving. 
Using  only  4  elements  and  a  small  number  of  nonlinear  iterations,  our  methods 
are  capturing  the  correct  deflections. 


Horizontal  End  Deflection  (m) 

Iteration  num. 

Fried  [27] 

Nonlinear 

Linearized 

1 

0.5000000 

0.46857929 

0.46857929 

2 

0.4337216 

0.41195049 

0.39038892 

3 

0.4124599 

0.41097754 

0.41736212 

4 

0.4109994 

0.41097726 

0.40889655 

5 

0.4109928 

0.41097726 

0.41164489 

6 

0.4109928 

0.41097726 

0.41076195 

Converged  Value 

0.4109928 

0.41097726 

0.41097726 

Table  4.1:  Horizontal  end  deflections  for  the  first  6  iterations  and  the  converged 
values  for  an  end  load  of  Fi  =  1.5.  The  results  shown  are  from  Fried  [27]  and 
the  nonlinear  hnite  element  method  and  linearized  method  presented  above. 

For  a  similar  problem  Fried  [27]  also  presents  convergence  data.  An 
end  load  Fi  =  5.0  is  used  with  the  number  of  elements  used  ranging  from  1 
to  9.  The  deflection  of  the  end  of  the  beam  is  the  value  of  interest.  Table 
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4.2  gives  the  end  deflection  for  Fried’s  method  and  our  nonlinear  method  for 
each  number  of  elements.  Our  linearized  method  failed  even  using  incremental 
loading  for  this  problem.  This  is  not  unexpected  for  such  a  large  deflection. 
For  each  refinement  level,  our  method  gives  a  smaller  relative  error  than  the 
previous  results.  To  calculate  the  relative  error  we  took  the  true  value  to  be 
0.713791524  m  which  is  the  deflection  from  our  nonlinear  method  using  100 
elements.  Fried’s  method  gives  a  convergence  rate  of  3.45  and  our  method 
gives  a  rate  of  4.13.  We  do  not  have  a  good  explanation  for  our  high  order 
of  convergence.  Using  a  7  element  configuration,  Fried  [27]  gives  deflections 
at  stable  equilibria  for  the  full  length  of  the  beam  for  end  loads  ranging  from 
0.5  to  9.0  in  increments  of  0.5.  We  performed  the  same  tests,  and  our  results 
can  be  seen  in  Figure  4.3.  Deflection  data  from  our  nonlinear  method  matches 
exceedingly  well  with  Figure  3  in  the  article  by  Fried  [27]. 

We  take  a  suite  of  benchmark  problems  for  distributed  loads  from  Dado 
and  Al-Sadder  [16].  For  plant  bending  due  to  fluid  resistance  distributed  loads 
are  a  much  bigger  factor  than  end  loads,  so  these  tests  are  incredibly  important 
to  verify  our  methods.  Dado  and  Al-Sadder  [16]  approximate  the  center  line 
of  the  beam  as  a  high-order  polynomial  and  formulate  a  method  to  minimize 
the  integral  of  the  residual  error.  They  verify  their  method  against  results 
from  the  MSC/NASTRAN  commercial  finite  element  package.  This  packages 
divides  the  beam  into  a  large  number  of  two-dimensional  beam  elements.  They 
perform  tests  for  prismatic  and  non-prismatic  beams,  but  since  our  interest  is 
only  for  prismatic  beams,  we  only  look  at  those  tests. 
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Figure  4.3:  Stable  deflection  configurations  for  a  cantilever  beam  loaded  with 
a  range  of  end  loads. 
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The  first  test  is  for  uniformly  distributed  loads  in  the  lateral  direction. 
The  configuration  for  this  test  is  shown  in  Figure  4.4.  A  range  of  uniform  dis¬ 
tributed  loads  q  =  —4,  —10,  —20,  —40,  —100  are  used.  Dado  and  Al-Sadder  [16] 
solved  for  the  stable  equilibria  using  their  method  with  a  15th  order  polynomial 
and  using  MSC/NASTRAN  with  100  elements.  MSC/NASTRAN  failed  due 
to  divergence  for  loads  of  —40  and  —100.  Figure  4.5  shows  stable  deflection 
conhgurations  given  the  designated  loads  using  our  fully  nonlinear  method.  10 
quadratic  beam  elements  and  a  nonlinear  tolerance  of  1.0  x  10“^  were  used. 
Essentially,  the  same  deflection  results  were  attained  using  the  nonlinear  hnite 
element  method  and  the  linearized  method  and  they  match  very  well  with  the 
results  from  Dado  and  Al-Sadder  [16].  There  was,  however,  significant  vari¬ 
ation  in  the  number  of  nonlinear  iterations  required  for  convergence  between 
the  two  methods.  The  number  of  iterations  also  depends  greatly  on  the  num¬ 
ber  of  incremental  loading  steps  and  the  plane  in  which  the  bending  occurs. 
Table  4.3  shows  the  total  number  of  Newton-Raphson  iterations  required  for 
each  load  value  using  the  fully  nonlinear  method  for  bending  in  the  x-z  plane, 
the  y-z  plane,  and  the  diagonal  plane.  We  see  that  using  more  than  1  load 
step  the  number  of  iterations  for  x-z  bending  and  y-z  bending  are  identical. 
However,  for  diagonal  bending  the  method  diverges  for  low  numbers  of  load 
steps.  We  also  see  divergence  for  very  large  deflection  in  the  y-z  plane  and 
only  one  step.  We  see  for  bending  in  each  direction,  that  taking  too  many 
incremental  loading  steps  can  cause  many  unnecessary  iterations  to  be  taken. 
Going  from  10  to  15  loading  steps  increases  the  total  number  of  iterations  for 
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bending  in  all  directions.  In  order  to  avoid  divergence  of  the  method  and  high 
computational  cost,  between  5  and  10  loading  steps  are  optimal. 

The  linearized  method  does  not  handle  large  deflections  nearly  as  well 
as  the  nonlinear  method.  Because  of  this,  we  could  only  test  this  method 
for  smaller  loads:  q  =  —2,  —4,  —6,  —8,  —10  Table  4.4  shows  the  total  number 
of  Picard  iterations  required  for  each  load  value  using  the  linearized  method. 
No  matter  how  many  loading  steps  were  used,  for  the  case  g  =  — 10  this 
method  required  over  1000  iterations  to  converge  for  bending  in  the  x-z  and 
y-z  planes.  It  failed  in  all  cases  for  q  <  —4  in  the  diagonal  plane.  This  ex¬ 
hibits  the  main  drawbacks  of  the  linearized  method.  It  works  relatively  well 
for  bending  in  planes  near  the  2  main  axes,  but  tends  to  fail  otherwise.  For 
very  small  deflections,  with  no  loading  it  converges  in  relatively  few  steps  and 
is  less  computationally  expensive  than  the  nonlinear  method.  However,  it  re¬ 
quires  a  prohibitive  number  of  Picard  iterations  for  larger  loads  or  diverges. 
Incremental  loading  does  not  seem  to  help  this  method,  as  in  each  case  do¬ 
ing  incremental  loading  increased  the  overall  cost.  While  we  have  seen  this 
method  work  relatively  well  for  end  loads,  it  is  not  particularly  successful  for 
distributed  loads. 

The  hnal  beam  benchmark  problem  from  Dado  and  Al-Sadder  is  one 
that  tests  if  the  methods  can  handle  truly  large  deflections.  The  conhguration 
is  shown  in  Figure  4.6  and  includes  a  uniform  vertical  load  q^  =  —5  and  a 
linearly  distributed  horizontal  load.  The  horizontal  load  g  is  0  at  the  base 
of  the  beam  and  qy  at  the  end.  We  take  values  of  g^  =  -l,-5,-10,-15,-20,- 
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Figure  4.4:  Configuration  of  the  uniformly  distrubted  load  test  problem. 


Deflections  Due  to  Uniform  Distributed  Loads 


Figure  4.5:  Deflections  for  the  uniformly  distrubted  load  test  problem. 
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Horizontal  End  Deflection  (m) 

Num.  of  Elements 

Eried  [27] 

Rel.  Error 

Nonlinear  Method 

Rel.  Error 

1 

0.7669329 

0.07444 

0.705610 

0.01146 

2 

0.718933 

0.00720 

0.713325 

0.00065 

3 

0.7151829 

0.001949 

0.713712 

0.00011 

4 

0.7143314 

0.00075 

0.7137673 

3.4E-5 

5 

0.7140374 

0.00034 

0.713782 

1.4E-5 

6 

0.7139174 

0.00017 

0.713787 

6.6E-6 

7 

0.7138622 

9.9E-5 

0.713789 

3.6E-6 

8 

0.7138340 

5.9E-5 

0.713790 

2.1E-6 

9 

0.7138185 

3.8E-5 

0.713791 

1.3E-6 

Conv.  Rate 

Table  4.2:  Convergence  results  for  an  end  load  Fi  =  5.0. 


Distributed  Loads 

-4 

-10 

-20 

-40 

-100 

bending  plane 

num.  of  load  steps 

Total  num.  of  iterations 

x-z  plane 

1 

4 

4 

33 

163 

100 

5 

15 

15 

17 

16 

40 

10 

30 

30 

30 

32 

31 

15 

42 

45 

45 

42 

43 

y-z  plane 

1 

4 

4 

30 

47 

* 

5 

15 

15 

17 

16 

34 

10 

30 

30 

30 

32 

31 

15 

42 

45 

45 

42 

43 

diagonal  plane 

1 

4 

5 

* 

* 

* 

5 

15 

15 

17 

17 

* 

10 

30 

30 

30 

32 

33 

15 

45 

45 

45 

47 

46 

Table  4.3:  Total  number  of  Newton-Raphson  iterations  for  the  uniform  dis¬ 
tributed  load  benchmark  problem  for  a  range  of  loads  in  different  planes  using 
the  fully  nonlinear  method.  *  denotes  that  the  method  diverged. 


30,-50,-80,  -100,  -120.  This  conhguration  leads  to  extremely  large  deflections. 
For  the  case  qy  =  —120  the  stable  equilibrium  has  the  beam  bent  nearly  180 
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Distributed  Loads 

-2 

-4 

-6 

-8 

-10 

bending  plane 

num.  of  load  steps 

Total  num.  of  iterations 

x-z  plane 

1 

4 

7 

13 

31 

1516 

5 

17 

24 

34 

59 

1483 

10 

30 

43 

60 

94 

1455 

20 

55 

77 

106 

162 

1497 

y-z  plane 

1 

4 

7 

13 

31 

1516 

5 

17 

24 

34 

59 

1483 

10 

30 

43 

60 

94 

1455 

20 

55 

77 

106 

162 

1497 

diagonal  plane 

1 

8 

29 

* 

* 

* 

5 

24 

46 

* 

* 

* 

10 

41 

71 

* 

* 

* 

20 

74 

119 

* 

* 

* 

Table  4.4:  Total  number  of  Picard  iterations  for  the  uniform  distributed  load 
benchmark  problem  for  a  range  of  loads  in  different  planes  using  the  linearized. 
*  denotes  that  the  method  diverged. 


degrees  from  the  unloaded  state.  Using  the  linearized  method,  these  very  large- 
deflection  problems  always  diverged.  Figure  4.7  shows  the  stable  equilibria 
positions  under  these  particular  loads  using  our  nonlinear  method.  These 
match  up  almost  perfectly  with  results  from  Dado  and  Al-Sadder’s  polynomial 
method.  The  MSC/NASTRAN  method  failed  to  converge  for  qy  <  —20,  but 
our  nonlinear  results  match  up  with  those  results  for  the  smaller  loads. 

Table  4.5  shows  the  total  number  of  Newton-Raphson  iterations  re¬ 
quired  for  convergence  for  bending  in  different  planes  and  with  different  num¬ 
bers  of  incremental  loading  steps.  We  see  identical  numbers  of  iterations  for 
each  case  in  the  x-z  and  y-z  planes.  It  is  interesting  that  for  these  cases,  even 
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Distributed  Loads 

-1  1  -5  1  -10  1  -15  1  -20  1  -30  1  -50  |  -80  |  -100  |  -120 

plane 

load  steps 

Total  num.  of  iterations 

x-z 

1 

4 

5 

22 

4* 

3* 

3* 

3* 

2* 

2* 

2* 

5 

15 

17 

19 

22 

26 

16* 

14* 

14* 

12* 

12* 

10 

29 

30 

35 

35 

37 

47 

26* 

45* 

25* 

to 

CO 

* 

25 

50 

69 

71 

75 

78 

81 

84 

* 

00 

F- 

96* 

78 

50 

100 

124 

136 

140 

142 

147 

146 

147 

147 

149* 

100 

200 

200 

255 

260 

256 

249 

249 

247 

248 

248 

y-z 

1 

4 

5 

22 

4* 

3* 

3* 

3* 

2* 

2* 

2* 

5 

15 

17 

19 

22 

26 

16* 

14* 

14* 

12* 

12* 

10 

29 

30 

35 

35 

37 

47 

26* 

45* 

25* 

23* 

25 

50 

69 

71 

75 

78 

81 

84 

■u 

00 

* 

96* 

78 

50 

100 

124 

136 

140 

142 

147 

146 

147 

147 

149* 

100 

200 

200 

255 

260 

256 

249 

249 

247 

248 

248 

diag. 

1 

4 

6 

* 

* 

* 

* 

* 

* 

* 

* 

5 

15 

16 

21 

* 

* 

* 

* 

* 

* 

* 

10 

30 

30 

36 

40 

* 

* 

* 

* 

* 

* 

25 

75 

75 

75 

83 

88 

102 

* 

* 

* 

* 

100 

200 

200 

259 

276 

282 

288 

317 

* 

* 

* 

200 

400 

400 

400 

510 

547 

567 

591 

621 

* 

* 

500 

1000 

1000 

1000 

1000 

1000 

1178 

1242 

1294 

1314 

1325 

Table  4.5:  Total  number  of  Newton-Raphson  iterations  for  mixed  distributed 
load  benchmark  problem  for  a  range  of  horizontal  linearly  distributed  loads.  A 
plain  *  denotes  that  the  method  diverged.  A  number  followed  by  *  denotes  that 
the  method  converged  to  an  unstable  equilibrium  in  that  number  of  iterations. 


for  the  very  large  loads  with  a  small  number  of  load  steps  the  method  always 
converged.  However,  as  seen  in  Table  4.5,  it  often  converged  to  an  unstable 
equilibrium.  These  are  not  states  that  are  global  minima  of  the  beam  energy 
functional,  so  they  are  not  the  correct  deflections.  Using  a  large  number  of  in¬ 
cremental  loads  allows  us  to  avoid  this  problem,  though.  With  100  load  steps 
the  method  converged  to  stable  equilibria  even  for  the  greatest  loads.  For  more 
reasonable  loads,  once  again  between  5  and  10  load  steps  were  computationally 
optimal.  For  bending  in  the  diagonal  plane,  rather  than  converging  to  unsta¬ 
ble  equilibria,  the  method  tended  to  diverge  if  too  few  incremental  load  steps 
were  taken.  The  method  also  required  many  more  total  iterations  to  converge 
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Figure  4.6:  Configuration  of  the  uniformly  distrubted  load  test  problem. 

in  the  diagonal  plane  than  the  other  planes.  Still,  taking  500  incremental  load 
steps,  it  converged  to  the  correct  stable  equilibria,  even  for  the  greatest  loads. 

4.3.5  Discussion  and  Conclusions  Concerning  Beam  Methods 

The  most  obvious  conclusion  from  the  benchmark  problems  it  is  that 
the  fully  nonlinear  method  of  minimizing  the  energy  functional  of  the  beam 
using  adequate  incremental  loading  is  an  extremely  robust,  computationally 
inexpensive,  highly  accurate  method  for  finding  the  stable  deflected  state  of  a 
cantilever  beam  with  an  end  or  distributed  load  applied  to  it.  This  method 
matched  results  from  several  other  methods  and  worked  even  for  incredibly 
large  deflections.  The  deflections  shown  in  4.7  are  much  greater  than  we  would 
expect  for  any  plant  interacting  with  a  realistic  fluid  flow.  For  reasonable  loads, 
generally  5  or  10  incremental  steps  were  enough  to  ensure  convergence  without 
taking  a  great  number  of  unnecessary  iterations.  We  have  also  shown  that 
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Deflections  Due  to  Mixed  Distributed  Loads 


Figure  4.7:  Deflections  for  the  uniformly  distrubted  load  test  problem. 

the  method  handles  complex  loads  in  more  than  one  direction.  The  method 
converged  with  many  fewer  iterations  if  the  deflection  occurs  close  to  the  x-z 
or  y-z  planes.  Away  from  those  planes,  taking  enough  incremental  load  steps 
ensures  convergence.  Therefore  for  flow  coupling  scenarios,  it  is  best  to  align 
one  of  those  planes  with  the  primary  flow  direction  if  known.  This  is  the 
direction  in  which  most  of  the  bending  will  occur,  so  it  should  cause  fewer 
iterative  steps  to  need  to  be  taken. 

The  method  converges  at  a  very  high  rate,  so  for  most  cases  only  a  small 
number  of  beam  elements  are  necessary.  For  modeling  flexible  vegetation  with 
this  method,  having  at  least  5  or  6  beam  elements  per  plant  and  using  10 
incremental  loading  steps  is  an  adequate  model  that  is  robust  and  unlikely 
to  diverge  for  a  realistic  load  due  to  flow  resistance.  When  modeling  a  larger 
scale  flow,  we  may  have  thousands  or  millions  of  plants  in  the  flow  modeled 
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as  cantilever  beams.  Therefore,  this  low  computational  cost  and  robustness  is 
essential.  Taking  too  many  iterations  and  divergence  issues  would  be  a  major 
bottleneck  in  a  coupled  flow/bending  code. 

The  linearized  method  is  much  more  problematic  than  the  nonlinear 
method.  We  see  that  it  handles  small  and  moderate  end  loads  very  well  without 
requiring  a  large  number  of  iterations.  It  is  important  to  note  that  each  Picard 
iteration  for  this  method  is  much  cheaper  than  a  Newton-Raphson  iteration 
for  the  previous  method.  However,  for  even  moderate  distributed  loads  this 
method  has  a  tendency  to  diverge  even  with  a  large  number  of  incremental  load 
steps  taken.  For  small  deflections  of  the  beam  the  linearized  method  is  just 
as  accurate  and  much  cheaper  computationally  than  the  nonlinear  method. 
However,  past  a  threshold  load  it  either  requires  a  huge  number  of  iterations 
to  converge  or  diverges.  It  is  important  to  note  that  the  “small”  deflections  for 
which  this  method  works  are  still  considered  to  be  “large”  deflections  in  beam 
modeling.  In  a  scenario  where  the  loads  are  expected  to  be  quite  small,  the 
linearized  method  is  cheaper  than  the  nonlinear  method  and  nearly  as  reliable. 
For  most  scenarios,  however,  the  nonlinear  method  is  much  more  reliable. 

4.4  Coupling  Flow  and  Flexible  Vegetation  Models 

4.4.1  Loads  and  Sink  Terms 

The  bending  of  the  vegetation  is  caused  by  a  distributed  load  that 
results  from  the  drag  force  caused  by  the  flow  moving  past  it.  The  concept 
of  drag  is  discussed  in  depth  in  Section  3.2.  The  distributed  load  q  is  a  force 
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per  unit  length  of  the  beam  in  units  kg/s^.  As  described  by  Kutija  and  Hong 
[53]  the  distributed  load  due  to  the  drag  of  a  fluid  with  density  p  and  mean 
velocity  V  =  [Vi,  V2,  V^]'^  over  the  length  of  the  beam  is 

qii-s)  =  Cdis)^b{s)pVi{s)^Vj{s)Vj{s)  (4.71) 

where  the  double  index  notation  implies  summation  in  the  standard  way,  Cd  is 
the  drag  coefficient  for  one  obstacle  and  b  is  the  diameter,  s  is  the  coordinate 
along  the  center-line  of  the  beam.  Hoerner  [36]  describes  drag  coefficients 
for  commonly  shaped  objects  over  a  large  ranges  of  Reynolds  numbers.  For  a 
single  upright  circular  cylinder,  ~  1.2  for  Re  between  1.0  x  10^  and  3.0  x  10^. 
Almost  every  flow  that  we  want  to  deal  with  lies  in  this  range  and  a  circular 
cylinder  is  a  good  approximation  of  the  shape  of  a  stem. 

The  drag  force  that  causes  the  bending  of  the  vegetation  must  be  con¬ 
sidered  in  the  flow  equations  as  well.  This  is  generally  done  by  including  it  as 
a  sink  term  in  the  momentum.  Kutija  and  Hong  [53]  show  that  the  correct 
sink  term  to  balance  the  forces  is 

/i(x)  =  -Cd^b{x.)p{x.)Vi{x.)^Vj{x.)Vj{^)  (4.72) 

for  one  stem.  They  propose  that  the  force  per  unit  unit  volume  for  a  packing 
of  stems  is 

F,(x)  =  NJ,  (4.73) 

where  Ny  is  the  number  of  stems  per  unit  horizontal  area  (in  m“^)  for  a  packing 
of  stems.  Li  and  Xie  [54]  adjust  this  by  using  drag  coefficients  that  factor  in 
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foliage  on  the  plant  and  the  shape  of  the  bent-over  plants.  However,  these 
studies  overlook  the  concept  of  bulk  drag  as  discussed  in  Section  4.1.  We 
propose  that  the  sink  term  for  a  packing  of  stems  is 

Fi{x.)  =  -Cd^6(x)Wp(x)l/i(x)y^Hj(x)l/j(x)  (4.74) 

where  Cd  is  not  necessarily  equal  to  Cd-  Note  that  this  force  is  equivalent  to  the 
beam  force  F  described  in  Section  4.2.1.  We  couple  our  fluid  and  structure 
models  using  these  drags  and  momentum  sinks  within  the  framework  of  an 
immersed  boundary  method. 

4.4.2  The  Immersed  Boundary  Method 

The  immersed  boundary  (IB)  method  is  a  method  for  modeling  the 
interaction  between  incompressible  Newtonian  fluids  and  elastic  solids.  It  was 
developed  by  Peskin  [77]  and  is  mostly  used  for  modeling  biological  systems, 
especially  cardiovascular  systems.  Generally,  fluid  models  are  defined  in  Eule- 
rian  frameworks  and  elastic  structure  models  are  dehned  in  Lagrangian  frame¬ 
works.  The  IB  method  allows  for  separate  Eulerian  and  Lagrangian  domains 
and  meshes.  Using  integral  transforms  involving  the  Dirac  delta  generalized 
function,  the  effects  of  the  Eulerian  domain  are  mapped  to  the  Lagrangian  do¬ 
main  and  vice  versa.  This  is  convenient  so  that  the  most  natural  frameworks 
may  be  used  for  the  fluid  and  structure  models,  and  simple  meshes  can  be  used. 
The  meshes  also  do  not  move,  as  they  do  in  the  other  most  common  method  of 
modeling  fluid-structure  interaction,  the  arbitrary  Lagrangian- Eulerian  (ALE) 
method.  In  ALE,  the  meshes  move,  and  for  large  deformations,  it  can  be 
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Figure  4.8:  The  separate  Lagrangian  and  Eulerian  domains  U  and  respec¬ 
tively  for  the  immersed  boundary  method,  x  maps  U  into  ^2. 

twisted.  The  IB  method  has  been  applied  to  flow  around  beams  by  Borges  et 
al.  [9];  however,  it  was  for  small-deflections. 

Suppose  U  is  the  Lagrangian  domain  for  the  elastic  structure  and  that 
n  is  the  Eulerian  fluid  domain.  In  the  Eulerian  sense,  U  changes  in  time.  Let 

X--iU,t)^n  (4.75) 

be  a  mapping  of  the  structure  domain  at  time  t  into  the  fluid  domain,  such 
that  x(s,t)  is  the  location  of  the  structure  at  coordinate  s  and  time  t  in  the 
fluid  domain  hi.  hi  \  x{U,t)  is  the  physical  region  occupied  by  the  fluid  and 
x{U,  t)  is  the  physical  region  occupied  by  the  structure  at  time  t.  Formulations 
for  the  IB  method  for  standard  linear  and  nonlinear  elasticity  problems  are 
described  thoroughly  by  Peskin  [77]  and  Griffith  and  Luo  [34]. 

The  fundamental  principle  behind  the  IB  method  is  that  one  may  use 
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integral  transforms  involving  the  Dirac  delta  to  map  between  the  Eulerian  and 
Lagrangian  frameworks.  Notice  that  (4.71)  and  (4.72)  involve  values  that  are 
naturally  functions  of  x  and  s.  When  modeling  flow  with  many  resolved  stems 
the  sink  term  (4.72)  should  not  be  imposed  everywhere  in  the  flow.  It  must 
only  be  imposed  where  there  is  a  stem.  Here  we  are  assuming  that  the  stem 
is  a  ID  manifold;  although,  the  3D  nature  of  it  is  taken  into  consideration  by 
the  stem  diameter  and  drag  coefficient  terms  in  (4.71)  and  (4.72).  We  suppose 
that  the  ID  beam  domain  is  immersed  in  the  3D  fluid  domain.  Using  the  IB 
method,  we  can  calculate  the  distributed  load  on  each  beam  and  momentum 
sink  per  unit  volume  on  the  fluid  by 

qi{s)  =  Cdh){s)  j  p(x)nj(x)^t;j(x)t;j(x)5(x  -  x(s))dx  (4.76) 

and 

Fi(x)  =  -C'rf^p(x)ni(x)^nj(x)nj(x)  j  b{s)5{l^  -  x{s))ds  (4.77) 

where  U  is  the  ID  center-line  of  the  the  beam  and  D  is  the  3D  fluid  domain. 
Note  that  these  formulas  involve  the  variable  velocity  v,  not  V,  the  mean 
over  the  length  of  the  beam.  By  integrating  over  the  length  of  the  beam  the 
mean  effect  is  captured,  x  maps  the  beam  location  in  terms  of  the  arc  length 
coordinate  s  to  the  location  in  the  fluid  domain  D.  x(s)  is  equivalent  to  r(s) 
introduced  in  Section  4.2  and  dehned  as  solutions  to  hrst  order  ordinary  differ¬ 
ential  equations  in  terms  of  Euler  angles  in  (4.29)  -  (4.31).  These  simple  ODEs 
may  easily  be  solved  by  numerical  integration  using  quadrature.  Therefore, 
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when  the  beam  system  is  solved  for  the  Euler  angles  given  a  certain  load,  x 
can  easily  be  calculated. 

The  above  formulation  ensures  that  the  momentum  for  the  coupled 
beam/fluid  system  is  conserved.  The  exact  force  that  goes  into  bending  the 
beam  is  lost  from  the  fluid  in  the  source  term.  We  see  this  by  integrating  the 
distributed  load  q  over  the  beam  and  integrating  the  sink  F  over  the  fluid 
domain: 


In  a  coupled  physical  fluid/beam  system  there  is  an  additional  load  on 
the  beams  due  to  weight  and  buoyancy.  Let  p/  be  the  density  of  the  fluid  and 
Ps  be  the  density  of  the  structure.  Assuming  a  beam  with  a  cross  sectional 
area  A,  the  load  on  the  beam  due  to  gravity  and  buoyancy  is 

=  {pf  -  Ps)A.  (4.79) 

In  the  beam  solving  step  gs,  must  be  added  to  q^;  however,  it  is  not  included 
in  the  above  force  balance  because  momentum  is  not  being  taken  away  from 
the  fluid  because  of  buoyancy. 
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4.4.3  Numerical  Concerns  for  Immersed  Bonndary  Method 


In  practice,  we  cannot  use  the  Dirac  delta  generalized  function  because 
we  integrate  using  quadrature  and  do  not  want  to  have  to  make  changes  to 
the  meshes.  The  quadrature  points  for  the  beams  and  the  fluid  would  have  to 
match  up  for  us  to  be  able  to  use  S,  which  would  be  quite  difficult  to  implement 
for  a  moving  immersed  structure.  Instead  we  use  an  approximation  to  the 
Dirac  delta  It  is  essential  that  this  approximation  is  bounded,  continuous, 
and  has  compact  support.  We  use  an  approximation  designed  for  the  immersed 
boundary  method  by  Peskin  [77]: 


,,(x)  = 


(4.80) 


where  h  is  the  beam  element  size  and 

fo, 


0(r)  =  < 


r  <  —2 

1(5  +  2r  -  V-7  -  12r  -  4r2),  -2  <  r  < -1 

-1  <  r  <  0 
0  <  r  <  1 


1(3  +  2r  +  Vl  -4r-4r2), 
|(3  —  2r  +  Vl  +  4r  —  dr^). 


(4.81) 


1(5  -  2r  -  V-7  +  12r  -  4r2),  l<r<2 

1,0,  2  <  r. 

The  function  0  is  shown  in  Figure  4.9.  For  this  formulation  of  6h,  we  see  that 

(5  as  h  — 0.  Also,  note  that  the  momentum  balance  (4.78)  still  holds 

replacing  S  with  Sh  and  integrating  over  both  B  and  D  using  quadrature.  Using 

6h  rather  than  6  leads  to  the  “smearing”  of  the  effects  of  the  Lagrangian  mesh 

on  the  Eulerian  mesh  and  vice  versa.  Since  we  are  interested  in  mean  flow 


characteristics,  this  is  not  a  problem,  since  momentum  is  conserved,  as  long 
as  the  meshes  are  properly  rehned. 
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Figure  4.9:  The  function  0  used  for  calculating  6h- 

We  see  that  using  approximation  to  6h  has  a  support  of  2  beam  element 
lengths  in  any  direction.  It  is  essential  that  the  mesh  spacing  is  chosen  in  such 
away  that  6h  is  a  good  approximation  to  6  for  integration  over  this  small 
region.  Griffith  and  Luo  [34]  investigate  mesh  spacing  for  integration  using 
Gaussian  quadrature.  If  h  is  the  mesh  spacing  for  the  Lagrangian  mesh  and 
hf  is  the  average  mesh  spacing  for  the  Eulerian  mesh,  they  experimented  with 
h/hf  =  1,  h/hf  =  2,  and  h/hf  =  4,  and  found  that  h/hf  =  2  gave  optimal 
results.  In  fact,  it  is  standard  practice  in  fluid-structure  interaction  problems 
to  use  a  Lagrangian  mesh  twice  as  coarse  as  the  Eulerian  mesh.  Therefore,  we 
choose  this  rehnement  ratio  for  all  simulations. 

We  verify  that  our  use  of  5h  leads  to  proper  calculations  of  loads  and 
sinks.  First  we  assume  the  fluid  has  a  uniform  velocity  prohle  Vi  =  1  m/s, 
V2  =  V3  =  0.  We  use  the  three  point  Gaussian  quadrature  rule  for  the  beam 
and  four  point,  second  order  Gaussian  quadrature  rule  for  the  3D  fluid  mesh. 
The  fluid  mesh  consists  of  unstructured  tetrahedra  with  a  maximum  element 
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diameter  half  that  of  the  uniform  element  size  of  the  ID  beams.  The  beam 
length  was  taken  to  be  0.5  m,  the  beam  radius  was  taken  to  be  0.01  m,  Cd  = 
1.2,  and  p  =  998.0  kg/m^.  Calculating  q  using  (4.76)  with  the  trne  5  and  then 
integrating  over  the  beam  should  give  a  total  drag  force  of  5.9889  kgm/s^. 
We  also  experimented  with  a  linear  velocity  prohle  Vi  =  z  with  all  other 
parameters  the  same.  For  this  case  the  total  drag  force  analytically  is  0.4990 
kgm/s^.  Table  4.6  shows  calculated  total  drag  values  and  the  relative  error 
for  flnid  mesh  sizes  ranging  from  0.1  m  to  0.0015625  m  for  both  nniform  and 
linear  velocity  prohles.  We  see  close  to  linear  convergence  rates  for  both  cases. 
For  mesh  rehnements  of  0.025  for  the  uniform  case  and  0.05  for  the  linear  case, 
we  see  error  less  than  5%  which  is  quite  reasonable  considering  the  complexity 
of  the  model. 


Vi  =  l 

Vi  =  z 

hf 

Total  Drag 

Relative  Error 

Total  Drag 

Relative  Error 

0.1 

5.0334 

0.159418 

0.5358 

0.07374 

0.05 

5.5497 

0.073196 

0.5201 

0.04228 

0.025 

5.7751 

0.035554 

0.5059 

0.01382 

0.0125 

5.8854 

0.017134 

0.5012 

0.00441 

0.00625 

5.9397 

0.008066 

0.4999 

0.001803 

0.003125 

5.9665 

0.003591 

0.4995 

0.001002 

0.0015625 

5.9799 

0.001353 

0.4995 

0.001002 

Analytical 

5.988 

0.4990 

Table  4.6:  Total  drag  force  valnes  calculated  using  immersed  boundary  formu¬ 
lation  for  nniform  and  linear  velocity  prohles. 
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4.4.4  Steady-State  Problem 

We  can  formulate  a  steady  fluid-structure  interaction  problem  using  the 
IB  formulation.  Rather  than  just  a  single  beam,  we  can  have  many  beams, 
say  ribeam-  Call  the  Lagrangian  domain  for  each  beam  and  let  bk  be  the 
diameter  function  and  Xk  be  the  mapping  into  hi  of  beam  k.  The  bulk  sink 
due  to  all  of  the  beams  is 

^beam 

f.(x)  =  5^  (4.82) 

k=l 

where 

= -C'rf^p(x)ni(x)y4^(x)n^y'  -  Xfc(s))ds.  (4.83) 

The  corresponding  distributed  load  on  beam  k  is 

(lfis)  =  Cd^bk{s)  j  p(x)ni(x)^nj(x)nj(x)(5,,(x  -  Xfc(s))dx.  (4.84) 

The  steady  flow  system  is 


V-v  =  0 

V  •  (v  (g)  V  —  z/(Vv Vv"^))  =  — -Vp-hg-h-F  (4.85) 

P  P 

with  the  sink  term  above  included.  We  can  develop  an  iterative  method  to 
solve  for  a  steady-state  solution  for  the  fluid-structure  interaction  problem; 

•  Initialize  fluid  system. 

•  Calculate  loads  on  beams  using  (4.84). 
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Figure  4.10:  Cantilever  beams  bent  to  steady-state  positions  due  to  uniform 
velocities  of  0.01  m/s  and  0.1  m/s  with  flexural  rigidity  El  =  3  x  10“®  kgm^/s^ 
solved  using  the  immersed  boundary  method. 

•  Solve  beam  systems  using  the  previous  step  as  initial  steps  for  iterative 
method.  Integrate  to  hnd  Xk- 

•  Calculate  sink  terms  using  (4.82). 

•  Solve  fluid  system  (4.85)  using  the  previous  step  as  initial  condition  for 
iterative  method. 

•  Repeat  until  both  fluid  and  beam  models  converge  within  given  tolerance 
from  step  to  step. 

This  gives  a  steady-state  solution  to  the  coupled  beam  and  fluid  model.  Figure 
4.10  shows  steady-state  bent  beam  prohles  caused  by  initial  uniform  velocity 
helds  of  0.01  m/s  and  0.1  m/s. 
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4.4.5  Dynamic  Problem 

For  the  dynamic  fluid-structure  interaction  problem,  we  utilize  a  quasi¬ 
static  model  for  the  beams  and  a  dynamic  Navier-Stokes  model  as  discussed 
in  Chapter  2.  Assume  that  we  have  time  steps  tn-  Now,  Xk)  F  are 

functions  of  time  as  well.  They  dynamic  flow  system  is 


V-v  =  0 

(9v  1  1 

^ V  ■  (v  0  V  —  z/(Vv Vv^))  =  — '\/p  +  g  +  -F.  (4.86) 

ot  P  P 

The  dynamic  fluid-structure  interaction  method  at  time  step  tn  is: 

•  Calculate  loads  on  beams  (4.84)  and  momentum  sinks  (4.82)  using  ve¬ 
locity  at  tn-l- 

•  Solve  beam  systems  using  the  previous  step  as  initial  steps  for  iterative 
method.  Integrate  to  And  Xki.'tn)- 

•  Solve  fluid  system  (4.86)  using  time  integration  to  time  tn- 

4.5  Model  Problems 

We  execute  the  immersed  structure  approach  for  two  types  of  model 
problems:  channel  flow  problems  and  wave  tank  problems.  These  two  types 
of  problems  were  chosen  for  a  variety  of  reasons.  Vegetated  channels  are 
extremely  common  in  nature.  A  vegetated  coastal  region  is  very  similar  to  a 
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vegetated  channel  with  an  extremely  large  horizontal  length.  Also,  there  have 
been  a  great  deal  of  experimental  studies  for  flow  through  vegetated  channels 
containing  rigid  and  flexible  vegetation.  We  can  reproduce  the  geometry  of 
these  channels  exactly  on  a  computational  domain  and  thus,  can  compare 
our  computational  results  to  experimental  results  for  the  exact  same  setup. 
We  compare  our  computational  results  to  experimental  studies  of  vegetated 
channels  by  Dunn  et  ah  [23]  and  Nepf  [73]  as  discussed  in  Section  4.1.2.  Also 
in  that  section,  we  discuss  the  connection  between  wave  attenuation  and  bulk 
drag.  The  existence  of  vegetation  has  a  major  effect  on  wave  attenuation  in 
near-shore  environments.  We  create  computational  “wave  tanks”  that  generate 
linear  planar  gravity  waves  using  a  two-phase  flow  method  and  measure  the 
mean  drag  as  they  travel  through  a  bed  of  vegetation.  We  reproduce  the 
geometry  from  experiments  by  Wu  et  al.  [115]  and  compare  calculated  bulk 
drag  coefficients  over  large  ranges  of  Reynolds  numbers.  It  is  important  to 
note  that  the  bulk  drag  coefficient  Cd  is  dependent  on  the  type  of  flow.  For 
a  channel  flow,  the  flow  travels  mostly  in  one  horizontal  direction,  passing 
through  the  bed  of  vegetation  as  it  flows.  For  a  wave-dominated  flow,  the 
fluid  moves  in  large  eddies  with  a  mean  velocity  in  a  horizontal  direction,  but 
passing  through  more  vegetation  as  it  rotates.  Therefore,  we  expect  the  bulk 
drag  coefficient  to  be  higher  for  a  wave  tank  than  for  channel  flow. 

As  in  these  experimental  studies,  we  arrange  the  vegetative  obstacles  in 
a  staggered  conhguration  as  shown  in  Figure  4.11.  The  distance  between  each 
stem’s  center  and  the  center  of  its  closest  neighbor  is  A.  The  corresponding 
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Figure  4.11:  Horizontal  spacing  of  vegetative  obstacles  for  channel  and  wave 
tank  setups. 

population  density  is  then  calculated  by 

Af„  =  Aa-2.  (4.87) 

We  assume  that  the  domains  have  smooth  bottoms,  so  that  there  is  no  bot¬ 
tom  friction.  To  model  rigid  vegetation,  we  use  a  variation  of  the  immersed 
structure  approach  discussed  in  the  previous  section.  Since  the  structures  will 
not  be  moving,  we  do  not  need  to  calculate  the  loads  on  them  or  recalculate 
the  mapping  Xk  i^^o  the  fluid  domain.  We  just  need  to  initialize  Xk  ^^^h 
beam  and  use  it  to  calculate  the  sink  term  at  each  fluid  quadrature  point  in 
the  fluid  model  at  each  time  step.  For  flexible  vegetation  we  use  the  full  dy¬ 
namic  immersed  boundary  fluid-structure  interaction  method  presented  in  the 
previous  section.  Beam  elasticities,  radii,  and  lengths  are  chosen  to  match  the 
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materials  used  in  the  experiments. 


Suppose  that  hi  is  the  vegetated  region  of  the  domain  in  which  we  are 
interested.  The  net  drag  force  per  unit  volume  is  then 

iFldx 

J  ’  4.88 

Jn  dx 

where  F  is  as  dehned  in  (4.1).  To  calculate  the  bulk  drag  coefficient  we  want 
the  mean  drag  force,  so  we  time  average: 

Fdrag  =  ,  ^  ,  [  Fdragdt,  (4.89) 

tf  -  ti  Jt. 

where  is  a  time  when  the  flow  is  first  developed  and  is  the  terminal  time. 
For  channel  flow,  Fdrag  tends  to  not  vary  a  great  deal  in  time,  so  this  arbitrary 
time  range  makes  little  difference.  We  calculate  the  corresponding  bulk  drag 
coefficient  by  using  the  bulk  drag  equation  (4.1): 


Cd  = 


2K 


drag 


pNjjbu"^ ' 


(4.90) 


4.5.1  Channel  Flow  Problem 

4. 5. 1.1  Problem  Setnp 

Our  model  channel  is  a  3D  rectangular  domain  4  m  long,  0.4  m  wide, 
and  0.368  m  tall.  We  drive  the  flow  by  imposing  an  inflow  velocity  Dirichlet 
boundary  condition  on  the  inflow  boundary.  The  inflow  velocity  is  positive  in 
the  x-direction  and  zero  for  the  y  and  z-directions.  For  velocity,  we  impose  a 
no-slip  condition  on  the  channel  bottom,  a  free  surface  condition  on  the  top  of 
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the  channel,  an  outflow  condition  on  the  outflow  boundary,  and  a  no-flux  con¬ 
dition  on  the  front  and  back  boundaries.  For  pressure,  we  impose  a  hydrostatic 
pressure  condition  on  the  outflow  boundary.  A  standard  acceleration  due  to 
gravity  g  =  9.8  m/s^  is  used.  The  water  density  used  is  p  =  998.2  kg/m^  and 
the  kinematic  viscosity  is  z/  =  1.004  x  10“®  m^/s  The  slope  of  the  domain  is 
0.  We  initialize  the  system  with  zero  velocity  and  hydrostatic  pressure  every¬ 
where  and  slowly  ramp  up  the  inflow  velocity  to  the  desired  value  u.  We  then 
run  the  model  out  in  time  using  second  order  backward  differentiation  formula 
(BDF)  time  integration  and  small  time  steps. 

The  model  vegetation  begins  1  m  from  the  inflow  boundary  of  the 
domain,  so  that  the  flow  is  able  to  fully  develop  before  entering  the  vegetated 
area.  The  vegetation  continues  to  the  outflow  boundary.  To  avoid  boundary 
effects  we  calculate  the  net  drag  only  in  the  center  of  the  channel,  the  region 
n  where 

n  =  [1.0, 4.0]  X  [0.1,  0.3]  X  [0,0.368].  (4.91) 

4.5. 1.2  Experimental  Studies 

Dunn  et  ah  [23]  conducted  their  experiments  in  a  19.5  m  long,  0.91  m 
wide,  0.61  m  deep  slanted  flume.  The  slope  of  the  flume  could  be  set  between 
0  and  10  percent.  Water  was  introduced  into  the  flume  from  a  constant  head 
tank  to  ensure  a  constant  discharge,  up  to  a  maximum  of  around  180  L/s. 
The  water  flowed  through  a  pattern  of  model  vegetation  in  a  pattern  similar 
to  that  in  Figure  4.11  and  exited  the  flume  into  a  large  sump  where  it  was 
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recirculated  into  the  flume.  For  the  experiments  that  we  are  comparing  to 
the  channel  slope  was  set  between  0.0036  and  0.0161.  Rigid  vegetation  was 
simulated  using  1/4  inch  diameter  wooden  dowels  6  inches  in  length.  Flexible 
vegetation  was  simulated  by  1/4  inch  diameter,  7  3/4  inch  long  plastic  drinking 
straws.  When  placed  onto  a  false  floor  the  average  height  of  the  rigid  dowels 
was  11.8  ±  1.67  cm  and  16.9  ±  1.61  cm  for  the  flexible  straws.  They  tested 
the  effect  on  the  bulk  drag  coefficient  with  respect  to  several  non-dimensional 
parameters:  Reynolds  number,  Fronde  number,  slope,  ratio  of  stem  diameter 
to  height,  and  ratio  of  stem  diameter  to  stem  spacing.  They  did  these  studies 
for  both  rigid  and  flexible  vegetation. 

Nepf  [73]  conducted  experiments  in  a  24  m  long  and  0.38  m  wide  glass- 
walled  flume.  The  flow  depth  H  =  0.15  m  was  always  used.  Vegetation 
was  modeled  by  using  0.064  m  diameter  wooden  dowels.  The  densities  of  the 
dowels  varied  between  200  to  2000  stems/m^.  This  is  the  range  of  actual  stem 
densities  observed  in  nature  by  saltwater  environments  by  Gambi  et  ah  [29] 
and  in  freshwater  environments  by  Kadlec  [43].  The  dowels  were  mounted  to 
Plexiglass  boards  and  randomly  placed  throughout  the  flume.  The  bulk  drag 
coefficient  was  calculated  over  for  various  population  densities  in  the  the  given 
range. 

4.5. 1.3  Results 

We  organized  our  channel  flow  computational  experiments  to  mimic 
the  geometry  of  the  experiments  by  Dunn  et  ah  [23].  For  rigid  and  flexible 
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experiments,  the  beam  length  is  0.1968  m,  the  beam  radius  is  0.003175  m, 
the  beam  flexibility  is  El  =  3.0  x  10“"^  Nm^  for  flexible  vegetation,  and  the 
population  density  Ny  is  172  stems/m^.  The  standard  mean  horizontal  velocity 
used  is  7/  =  0.57  m/s.  The  parameters  listed  here  are  the  standard  parameters. 
We  ran  several  simulations  holding  all  parameters  but  one  constant  and  varying 
one  to  study  the  effects  of  that  parameter  on  the  bulk  drag  coefficient. 

We  see  how  the  bulk  drag  coefficient  varies  with  respect  to  four  dimen¬ 
sionless  quantities  introduced  by  Dunn  et  ah  [23].  The  first  quantity  is  the 
Reynolds  number  Re  defined  with  the  channel  height  H  as  the  length  param¬ 
eter:  Re  =  UH/v.  This  is  a  dimensionless  measure  of  velocity.  The  second 
quantity  is  Dy/H  where  Dy  is  the  stem  diameter.  Here  we  keep  H  constant 
and  vary  Dy.  This  is  a  dimensionless  measure  of  obstacle  thickness.  The  third 
quantity  is  HDy/ \^.  We  keep  H  and  Dy  constant  and  vary  the  spacing  A.  This 
is  a  dimensionless  measure  of  population  density.  The  final  quantity  is  hy/H 
where  we  keep  H  constant  and  vary  hy.  This  is  a  non-dimensional  measure  of 
stem  height. 

Dunn  et  ah  [23]  introduce  2  versions  of  the  bulk  drag  coefficient.  The 
two  versions  differ  depending  on  the  vertical  length  over  which  the  bulk  drag 
coefficient  is  calculated.  For  CdH  the  vertical  scale  used  in  D  is  77  the  channel 
depth.  For  CdK  the  vertical  scale  used  in  D  is  hy  the  stem  height.  Normally, 
we  are  interested  mostly  in  Cdh^)  but  to  compare  to  their  results  we  calculate 
both. 

Figure  4.12  shows  how  the  the  four  dimensionless  parameters  effect 
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the  bulk  drag  coefficients  for  rigid  and  flexible  vegetation.  The  bulk  drag 
coefficients  calculated  by  Dunn  et  ah  [23]  are  much  more  scattered  than  ours. 
This  difference  is  most  likely  due  to  error  caused  by  them  having  to  measure 
velocity  information  experimentally  and  calculate  the  drag  coefficient  with  a 
model,  whereas  with  our  model,  we  can  calculate  the  net  drag  force  directly. 
While  our  results  are  not  as  scattered  as  theirs,  the  bulk  drag  coefficients  are 
in  the  same  ranges  and  have  similar  tendencies.  The  experimental  results 
show  CdK  varying  between  0.8  and  1.5  in  all  cases  and  usually  being  near 
1.0.  We  see  our  calculations  of  Cdh^  almost  always  being  between  0.9  and  1.0. 
This  matches  up  extremely  well.  The  experimental  calculations  of  Cdu  vary 
between  0.05  and  0.6,  mostly  being  between  0.1  and  0.4.  Our  calculations  are 
mostly  near  0.5.  Our  method  tends  to  calculate  only  slightly  higher  values  of 
CdH  than  the  experimental  results. 

The  most  noticeable  pattern  is  that  in  all  cases,  the  flexible  vegetation 
resulted  in  slightly  lower  bulk  drag  coefficients.  The  experimental  results  also 
show  this  basic  tendency.  Also,  the  bulk  drag  coefficients  tend  to  decrease 
slightly  with  increasing  Reynolds  number.  Bent  stem  prohles  for  a  range  of 
Reynolds  numbers  are  displayed  in  Figure  4.13.  Our  results  also  show  the 
bulk  drag  coefficient  being  mostly  constant  with  respect  to  changing  obstacle 
diameter  and  height.  Dunn  et  al.  [23]  also  came  to  that  conclusion.  The  steep 
slope  of  CdH  due  to  changing  length  is  due  to  how  we  calculate  this  value, 
and  it  does  not  indicate  an  increase  in  drag  due  to  increased  stem  length. 
Finally,  we  see  a  general  tendency  for  the  bulk  drag  coefficient  to  decrease 
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as  the  population  density  is  increased.  Dunn  et  al.  [23]  also  came  to  the 
conclusion  that  the  bulk  drag  coefficient  decreases  with  density.  They  found 
a  slightly  more  pronounced  relationship  than  our  results  show. 

These  results  show  that  our  immersed  structure  method  for  modeling 
flow  through  vegetated  regions  does  a  reasonable  job  at  calculating  bulk  drag 
coefficients.  Quantitatively,  we  see  that  we  calculate  CdK  ^  range  that 
matches  up  extremely  well  with  experimental  results.  Our  calculations  of 
CdH  are  only  slightly  higher  than  those  from  experiments.  Qualitatively,  our 
method  predicts  tendencies  of  the  bulk  drag  coefficient  to  vary  with  respect 
to  different  parameters  that  was  also  seen  by  experiments.  We  show  that  the 
drag  coefficients  tend  to  decrease  with  increased  velocity  and  with  increased 
population  density.  Rigid  obstacles  tend  to  have  higher  bulk  drag  coefficients 
than  flexible  ones.  Stem  thickness  and  length  do  not  have  major  impacts  on 
the  bulk  drag  coefficient.  That  our  model  captures  these  tendencies  that  have 
been  observed  in  experimental  settings  serves  to  vindicate  this  model. 

Nepf  [73]  only  studies  the  relationship  between  population  density  and 
the  bulk  drag  coefficient  scaled  with  the  vegetation  height,  Ch^.  As  with 
the  previous  experiments,  we  reproduce  the  geometry  of  the  staggered  array 
of  rigid  dowels.  We  ran  computational  experiments  with  stem  population 
densities  ranging  from  200  to  2000  stems/m^,  the  range  seen  in  nature,  as 
mentioned  before.  The  dimensionless  quantity  for  population  used  for  these 
experiments  is  as  in  Nepf.  Varying  this  dimensionless  quantity  from  0.003 
to  0.3,  we  chose  10  values  on  a  logarithmic  scale.  A  mean  inflow  velocity  oiU  = 
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Figure  4.12:  The  effect  of  dimensionless  values  representing  velocity,  stem  di¬ 
ameter,  stem  spacing,  and  stem  length  on  the  bulk  drag  coefficient  for  channel 
flow. 
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Figure  4.13:  Bent  vegetation  profiles  for  channel  flow  of  Reynolds  numbers  (a) 
0,  (b)  6.00  X  10^  (c)  7.17  x  10^  (d)  8.58  x  10^  (e)  1.03  x  10^  (f)  1.23  x  10^ 
(g)  1.47  X  10^,  and  (h)  1.75  x  10^  respectively. 
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0.57  m/s  was  used  in  all  cases.  Figure  4.14  shows  the  relationship  calculated 
between  dimensionless  population  density  and  the  bulk  drag  coefficient.  These 
results  match  up  exceptionally  well  with  the  experimental  results.  The  range 
of  population  densities  simulated  here  is  much  larger  than  those  done  when 
reproducing  the  experiments  of  Dunn  et  al.  [23].  The  decrease  of  the  bulk 
drag  coefficient  with  increased  population  density  is  much  more  pronounced 
over  this  range.  As  with  the  results  from  Nepf,  for  dimensionless  population 
density,  between  0.003  and  0.02  the  bulk  drag  coefficient  calculated  is  nearly 
constant  near  a  value  of  1.0.  Higher  than  this  population  density,  Ch^  drops 
in  a  gradual  arc  with  a  shape  very  similar  to  the  one  from  the  experiments. 
For  the  highest  population  densities,  Ch^  reaches  a  value  near  0.5.  Nepf  [73] 
saw  a  slightly  lower  value,  closer  to  0.4.  The  experimental  results  due  to 
Nepf  are  highly  cited  and  very  well-respected.  The  results  from  that  study 
are  compared  to  and  match  up  well  with  experimental  results  due  to  Seginer 
et  al.  [94],  Kays  and  London  [45],  Petryk  [79],  and  Zdravkovich  [117].  The 
fact  that  the  curve  of  bulk  drag  coefficient  versus  population  density  using 
our  immersed  structure  method  matches  the  curve  from  experimental  results 
extremely  closely  does  a  great  deal  to  verify  our  model.  It  shows  that  the 
immersed  structure  method  captures  the  fact  that  the  bulk  effect  of  a  dense 
packing  of  N  vegetative  obstacles  is  not  equivalent  to  N  times  the  effect  of 
one  vegetative  obstacle.  The  more  densely  the  stems  are  packed  past  a  certain 
threshold,  each  stem  has  less  of  a  drag  effect  on  the  flow. 
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Figure  4.14:  The  relationship  between  population  density  of  stems  and  the 
bulk  drag  coefficient. 

4.5.2  Wave  Tank  Problem 
4.5.2. 1  Two-Phase  Flow  Method 

To  implement  a  wave  tank  problem,  we  used  a  two-phase  flow  method. 
We  use  a  conservative  level  set  method  presented  by  Kees  et  ah  [46].  This 
method  is  designed  for  free-surface  flow  computations  with  complex  phenom¬ 
ena  like  waves.  It  uses  standard  level  set  methods  with  additional  discrete 
conservation  properties.  The  flow  model  is  the  same  LES  model  as  discussed 
in  Section  2;  however,  now  instead  of  one  incompressible  fluid  with  density  p 
and  kinematic  viscosity  z/,  we  have  a  water  phase  with  density  pw  and  viscosity 
Vw  and  an  air  phase  with  density  pa  and  Va-  We  take  pa  =  1.205  kg/m^  and 
Va  =  1.5  X  10“®m^/s.  The  water  parameters  are  the  same  as  in  the  channel 
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problem.  The  boundary  between  the  air  and  water  phase,  T  is  represented  by 
the  zero  level  set  of  a  function  0,  i.e. 

r  =  {x|0(x  =  0}.  (4.92) 

As  with  standard  level  set  methods,  the  interface  evolution  is  described  by  the 
level  set  equation: 

^+v-V0  =  O,  (4.93) 

where  v  is  the  fluid  velocity.  The  water  lies  in  the  region 

riu,  =  {x|0(x)  <  0},  (4.94) 

and  the  air  lies  in  the  region 

ila  =  {x|0(x)  >  0}.  (4.95) 

Using  the  method  formulated  by  Kees  et  ah  [46],  (4.93)  is  coupled  with  a 
volume  fraction  equation,  a  redistancing  equation,  and  a  level  set  mass  con¬ 
servation  equation  creating  a  conservative  method.  The  resdistancing  step 
solves  an  eikonal  equation  which  makes  0  be  a  signed  distance  function.  The 
density  and  viscosity  near  the  interface  are  made  continuous  by  using  a  regu¬ 
larized  Heaviside  function 

{0  (j)<-e 

Pl  +  ^  +  lsmy)  |.#.|<€  (4.96) 

1  4>>e 

where  e  is  a  small  parameter.  We  take  e  to  be  1.5  times  the  fluid  element 
diameter.  Thus,  p  and  v  are  dehned  by 

p  =  p^H,{<t>)+pa{l-H,m,  (4.97) 
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and 

u  =  iy^H,{(j))  +  Uail  -  H,{(j))).  (4.98) 

It  is  important  to  note  that  p  as  defined  above  is  the  one  used  in  the  drag 
calculations. 


4. 5. 2. 2  Linear  Plane  Wave  Theory 


The  wave  tank  experiments  by  Wu  et  al.  [115]  that  we  reproduce 
computationally  use  a  paddle  wavemaker.  Let  S  be  the  stroke  of  the  paddle, 
T  be  the  period  of  the  wave  being  generated,  L  be  the  wavelength  of  the 
wave  generated,  and  k  =  27^ j L  be  the  wavenumber.  Suppose  the  height  of  the 
undisturbed  wave  tank  is  h  and  the  coordinate  is  zero  at  this  height.  Suppose 
the  wave  amplitude  is  A.  For  progressive  waves,  the  wavemaker  frequency  a 
is  defined  by 

=  —gk  tan  kh  (4.99) 


where  g  =  —9.8  m/s^  is  the  gravitational  constant.  The  wave  period  is  T  = 
27r/(T.  Let  rj  be  the  wave  height  as  a  function  of  position  and  time  and  ^  be 
the  velocity  function.  The  governing  equation  for  a  linear  plane  gravity  wave 
in  the  x-z  plane  is  the  Laplace  equation 


dx‘^ 

where  ip  is  the  velocity  potential, 
boundary  conditions  are 


d^tP 


(4.100) 


The  dynamic  and  kinematic  free  surface 


V 


1  dip 
9  dt 


z  =  0, 


(4.101) 
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dip  df] 

respectively.  The  bottom  boundary  condition  is  a  no-flow  condition: 


(4.102) 


n  h 

=0  z  =  -h. 

dz 


(4.103) 


The  linearized  lateral  boundary  condition  is 


4i(0,i:,t)  =  ^  a  cos  at. 


(4.104) 


Using  standard  plane  wavemaker  theory  as  presented  by  Dean  and  Dalrymple 
[20],  the  analytical  solution  to  the  Laplace  equation  with  the  given  boundary 


conditions  gives 


rj{x,  t)  =  A  cos{kx  —  at), 


^  ,  —qkAcosh(k(z  +  h))  .  ,  ,  , 

^lix,  z,t)  = - sin  (at  -  /ex), 

a  cosh{kh) 

^2{x,Z,t)  =  0, 

^  ,  —qkAcosh(k(z  +  h))  ,  ,  , 

6(a^,  i:,  t)  = - cos{at  -  kx). 

a  cosh{kh) 


(4.105) 

(4.106) 

(4.107) 

(4.108) 


4. 5. 2. 3  Model  Setup 

We  use  the  above  plane  wave  theory  to  set  the  initial  conditions  and 
boundary  conditions  for  our  two-phase  flow  method.  The  wave  tank  domain 
that  we  use  is  shown  in  Figure  4.15.  The  tank  consists  of  a  rectangular  section 
with  a  wedge  section  on  the  end.  Let  the  rectangular  section  have  the  dimen¬ 
sions  Li  X  L2  X  L^.  The  wedge  section  is  an  absorbing  region  that  completely 
dampens  out  the  wave  to  prevent  reflections  off  the  right  boundary.  Given  a 
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Figure  4.15:  Domain  for  the  wave  tank  problem. 
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wave  amplitude  A,  wavelength  L,  and  initial  mean  flow  height  h,  the  level  set 
function  0  has  the  initial  condition 

(j)o{x,y,  z)  =  z  —  h  +  ri{x,0)  (4.109) 

the  velocity  initial  conditions  are  set  by 

vo{x,y,z)  =  ^{x,z  -  h,0)(l  -  (4.110) 

and  the  pressure  initial  condition  is  hydrostatic: 

Po{.x,y,z)  =  pgz.  (4.111) 

The  continued  propagation  of  this  plane  gravity  wave  is  forced  by  a  Dirichlet 
boundary  condition  on  the  left  inflow  boundary  of  the  wave  tank  domain.  This 
condition  is 

v(0,|/,^,t)  =  i{x,z-  h,t)(l  -  H,{(t))).  (4.112) 

A  hydrostatic  pressure  Dirichlet  condition  is  imposed  at  the  top  of  the  domain. 
A  no-flux  condition  is  imposed  on  all  other  boundaries.  In  the  wedge  on  the 
far  right  of  the  wave  tank  domain,  a  large  momentum  sink  term  is  imposed. 
This  term  is  large  enough  to  almost  completely  dampen  out  the  wave  before 
hitting  the  right  boundary.  This  minimized  the  reflection  of  the  waves  back 
through  the  tank.  This  is  designed  to  mimic  the  wave  absorbers  used  at  the 
end  of  the  tanks  used  by  Wu  et  ah  [115]. 

The  wave  tank  domain  used  for  simulations  had  length  Li  =  4.5  m, 
width  L2  =  0.69  m,  and  height  L3  =  1  m.  The  vegetated  area  began  at 
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X  =  0.5  m,  and  the  wave  absorbing  wedge  began  at  a;  =  3.5  m.  The  area 
between  these  lengths  is  hlled  with  vegetative  obstacles.  All  wave  simulations 
were  done  holding  the  wavelength  L  =  2.0  m  constant.  The  experiments 
done  by  Wu  et  ah  [115],  use  variable  wavelengths,  but  the  values  they  used 
are  all  close  to  this.  We  ran  many  simulations  with  each  type  of  vegetation 
for  different  wave  amplitudes  and  several  mean  heights.  In  order  to  get  an 
accurate  mean  bulk  drag  force  calculation,  we  run  several  waves  through  the 
vegetated  bed.  When  calculating  Fdrag  with  (4.89),  it  is  essential  that  tf  and 
ti  occur  at  the  same  place  in  the  wave  period,  i.e.  is  an  integer. 

4. 5. 2. 4  Experimental  Setup 

As  mentioned  numerous  times,  our  computational  wave  tank  setup  is 
designed  to  model  an  experimental  study  by  Wu  et  al.  [115].  Their  exper¬ 
imental  setup  included  a  20.6  m  long,  0.69  m  wide,  1.22  m  tall  wave  flume. 
A  paddle  wave  generator  created  regular  and  irregular  gravity  waves  on  one 
end,  and  a  porous,  parabolic  wave  absorber  was  placed  on  the  opposite  end  to 
minimize  reflections.  They  performed  studies  with  flat  and  slanted  bottoms. 
We  are  only  concerned  with  the  flat-bottom  experiments  for  now.  Between 
the  lengths  of  12.9  m  and  16.56  m  arrays  of  model  vegetation  were  arranged 
in  a  pattern  exactly  like  that  in  Figure  4.11.  They  utilized  several  types  of 
model  “vegetation”  of  which  we  focus  on  four.  The  first  is  a  rigid  vegetation 
type  using  wooden  dowels  similar  to  the  channel  flow  studies  by  Dunn  et  al. 
[23]  and  Nepf  [73].  These  dowels  have  a  diameter  of  9.4  mm  and  are  0.48  m 
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tall.  They  are  spaced  with  A  =  57.4  mm,  with  the  corresponding  population 
density  =  350  stems/m^.  Second,  their  flexible  “vegetation”  consists  of 
Ethylene  propylene  diene  monomer  (EPDM)  foam  rubber  cylinders,  also  with 
a  diameter  of  9.4  mm  and  height  of  0.48  m.  The  spacing  and  density  are  the 
same  as  the  rigid  model  as  well.  These  foam  hoses  have  a  density  of  368  kg/m^ 
and  Young’s  modulus  E  =  4  MPa. 

The  last  two  types  of  model  vegetation  are  real  beds  of  Spartina  al- 
terniflom,  an  extremely  common  type  of  coastal  marsh  grass.  One  of  the  beds 
of  grass  is  dormant  and  one  is  green.  The  plants  were  grown  in  an  outdoor 
nursery  and  transported  to  the  experimental  facility.  The  stem  density  was 
determined  to  be  Ny  =  545  stems/m^  for  the  dormant  grass  and  Ny  =  405 
stems/m^  for  the  live  grass. 

We  model  the  rigid  vegetation  in  the  same  manner  as  in  the  channel 
flow,  using  the  immersed  structure  method  without  allowing  bending.  The 
foam  rubber  hose  is  modeled  as  flexible  cantilever  beams  with  the  fully  non¬ 
linear  beam  solving  method.  Young’s  modulus,  E,  is  known  and  the  second 
moment  of  area  is  that  of  a  circular  cylinder  I  =  =  7.66  x  10“^°  We 

space  the  stems  in  the  same  manner  as  with  the  channel  flow,  using  the  value 
of  A  shown  in  Figure  4.11. 

To  model  the  real  vegetation,  we  also  model  each  plant  as  a  flexible 
cantilever  beam.  For  the  parameters  describing  the  stem  we  use  experimental 
data  from  Wu  et  al.  [115].  They  found  that  for  dormant  and  green  grass,  a 
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good  approximation  for  the  modulus  of  elasticity  is 

/  ,  N  1.773 

E  =  183204  l^jyj  ,  (4.113) 

where  is  the  vegetation  diameter,  is  the  vegetation  height,  and  the  units 
for  E  are  N/m^.  Different  values  were  found  by  Hessini  et  ah  [35].  Wu  et 
al.  [115]  also  measured  other  parameters.  Samples  of  Spartina  alterniflora 
were  taken  from  held  locations  in  Louisiana  and  Mississippi,  and  stem  height, 
length,  density,  and  plant  volume  were  measured  for  each  stem.  They  give 
mean  values  and  standard  deviations  for  stem  diameter  and  height.  We  calcu¬ 
late  a  truncated  normal  distribution  based  on  the  mean,  standard  deviation, 
maximum,  and  minimum  values.  We  chose  a  truncated  normal  distribution 
rather  than  a  normal  distribution  because  the  outliers  of  the  normal  distri¬ 
bution  were  physically  unrealistic.  Figure  4.16  shows  the  probability  density 
functions  for  stem  length  and  radius  for  dormant  Spartina  alterniflora,  and 
Figure  4.17  gives  the  same  information  for  green  Spartina  alterniflora.  The 
mass  density  for  all  Spartina  alterniflora  was  close  to  610  kg/m^  so  that  value 
was  chosen  for  all  stems. 

Using  (4.87),  and  the  known  population  densities  iV„  for  the  real  veg¬ 
etation,  we  calculate  the  spacing  A.  We  arrange  the  plants  in  our  standard 
array  pattern  with  this  spacing.  For  each  stem  in  the  array,  we  choose  a  radius 
and  length  from  the  probability  distributions  presented  above.  We  use  the  ra¬ 
dius  to  calculate  the  second  moment  of  area  /.  Hence,  this  provides  us  with 
a  bed  of  vegetation  that  is  somewhat  representative  of  a  real  bed  of  Spartina 
alterniflora. 


138 


159 


Figure  4.16:  Length  and  radius  probability  density  functions  for  dormant 
Spartina  alterniflora. 


Figure  4.17:  Length  and  radius  probability  density  functions  for  green  Spartina 
alterniflora. 
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4. 5. 2. 5  Results 

As  with  the  channel  flow  problems  we  only  calculate  the  drag  force 
Fdrag  over  the  center  of  the  domain  to  avoid  boundary  effects.  We  used  a 
highly-refined  unstructured  tetrahedral  mesh.  The  maximum  element  diame¬ 
ter  allowed  is  0.025  m.  The  resulting  meshes  contain  approximately  2  million 
elements.  We  run  at  least  5  waves  through  the  vegetated  part  of  the  domain 
and  time  average  to  calculate  the  mean  drag  force  Fdrag  and  then  the  bulk 
drag  coefficient  Cd-  For  these  simulations,  the  vertical  length  scale  for  the 
region  over  which  the  calculation  of  the  bulk  drag  is  performed  is  the  mean 
undeflected  vegetation  height  hy.  Therefore  it  is  similar  to  CdK  presented  in 
the  channel  flow  section. 


The  average  velocity  u  used  in  calculating  the  bulk  drag  coefficient 
with  (4.1)  is  the  time-averaged,  depth-averaged,  magnitude  of  the  horizontal 
velocity.  This  is  calculated  analytically  as 

fi  =  \a^,  (4.114) 


where  A  is  the  wave  amplitude  and  h  is  the  mean  water  depth.  The  Reynolds 
number  Re  is  defined  using  the  mean  stem  diameter  Dy  as  a  length  scale  and 
the  maximal  horizontal  velocity  just  before  the  vegetation  zone  at  mid-height 


of  the  plants 


Ur  =  Aa 


cosh.{ka*h) 


(4.115) 


sinh(A;a*h)’  '  ‘  ' 

where  a*  is  hy/h  for  submerged  vegetation  and  unity  for  emergent  vegetation. 
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Therefore, 

Re  =  (4.116) 

For  each  type  of  model  vegetation,  simulations  were  run  for  three  dif¬ 
ferent  mean  water  depths.  For  the  rigid  and  model  flexible  vegetation,  these 
depths  were  0.4  m,  0.5  m,  and  0.6  m.  This  includes  one  emergent  vegetation 
example,  one  submerged  vegetation  example,  and  one  example  where  the  wa¬ 
ter  depth  is  nearly  the  unbent  vegetation  height.  For  the  cases  modeling  real 
vegetation  depths  of  0.5  m,  0.6  m,  and  0.7  m  were  used,  also  giving  one  emer¬ 
gent  vegetation  example,  one  submerged  vegetation  example,  and  one  example 
where  the  water  depth  is  nearly  the  unbent  vegetation  height.  We  plot  calcu¬ 
lated  bulk  drag  coefficients  using  our  immersed  structure  method  for  specihc 
Reynolds  numbers  for  the  3  different  depths  and  a  curve  £t  for  the  results  due 
to  Wu  et  ah  [115]. 

Figure  4.18  shows  the  results  for  the  rigid  vegetation  model.  The  curve 
fit  to  the  SERRI  data  was  calculated  in  the  report  as 

Cd  =  2.557  +  {565/ Ref 

We  see  that  for  very  low  Reynolds  numbers,  we  underestimate  Cd-  This  is  most 
likely  because  we  are  using  Cd  =  1-2  for  the  drag  coefficient  of  a  single  beam  in 
the  immersed  structure  model.  For  Reynolds  numbers  this  low,  this  is  probably 
a  smaller  value  of  Cd  than  what  is  realistic  in  nature,  causing  our  method  to 
result  in  a  lower  value  for  the  drag  coefficient.  Around  a  Reynolds  number  of 
500  and  higher,  our  results  match  up  very  closely  to  the  experimental  results. 
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Figure  4.18:  Bulk  drag  coefficient  values  for  rigid  vegetation  setup  compared 
to  Wu  [115]  results. 

We  capture  the  gradual  decrease  with  respect  to  Re  that  they  observed,  with 
a  similarly  shaped  dip.  As  with  the  the  experimental  results,  we  see  Cd  begin 
to  level  off  near  Reynolds  number  1500.  Past  that  we  calculate  slightly  lower 
values  than  them;  although,  their  is  a  large  spread  to  their  data  in  that  range. 
It  appears  that  from  these  results  that  mean  flow  depth  does  not  have  a  major 
impact  on  the  bulk  drag  coefficient. 

Figure  4.19  shows  a  time  series  of  a  wave  traveling  through  a  bed  of 
model  flexible  vegetation.  Notice  how  our  immersed  structure  method  cap¬ 
tures  the  complex  bending  back  and  forth  that  occurs  as  a  water  wave  travels 
through  the  domain.  Figure  4.20  shows  the  bulk  drag  coefficient  results  for 
the  foam  rubber  model  flexible  vegetation  model.  The  curve  £t  to  the  SERRI 
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data  was  calculated  in  the  report  as 

Cd  =  (4020/i?e)°-®®h 

As  with  the  rigid  model,  we  underestimate  the  bulk  drag  coefficient  greatly  for 
very  small  Reynolds  numbers.  The  explanation  for  this  is  probably  the  same 
as  for  the  rigid  vegetation.  Around  a  Reynolds  number  of  400,  our  calculations 
begin  to  match  up  exceedingly  well  with  the  experimental  results.  We  capture 
the  decreasing  slope  with  increasing  Reynolds  number.  Our  data  points  stay 
almost  exactly  on  the  best  £t  curve  for  the  experimental  results  up  until  a 
Reynolds  number  of  nearly  2500.  That  our  calculations  match  up  so  closely 
to  experimental  results  for  such  a  complicated  problem  is  incredibly  exciting. 
This  does  much  to  verify  our  proposed  method  of  modeling  flexible  vegetation 
for  calculating  the  bulk  drag  coefficient.  Once  again,  it  appears  that  from 
these  results  that  mean  flow  depth  does  not  have  a  major  impact  on  the  bulk 
drag  coefficient  for  flexible  vegetation. 

Our  results  for  the  models  approximating  real  Spartina  alterniflora 
stems  were  not  nearly  as  good  as  those  for  the  rigid  and  flexible  model  vegeta¬ 
tion.  For  both  the  model  dormant  and  green  grass,  we  greatly  underestimated 
Cd  for  low  Reynolds  numbers  as  with  the  other  models.  The  model  dormant 
Spartina  alterniflora  data  matches  the  experimental  data  relatively  well  start¬ 
ing  near  a  Reynolds  number  of  500.  The  curve  fit  to  the  SERRI  data  was 
calculated  in  the  report  as 

Cd  =  {2623/ 
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Figure  4.19:  A  wave  traveling  through  a  bed  of  model  flexible  vegetation. 
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Figure  4.20:  Bulk  drag  coefficient  values  for  flexible  vegetation  setup  compared 
to  Wu  [115]  results. 

From  our  data,  we  see  Cd  leveling  off  here  and  becoming  almost  constant. 
The  experimental  data  shows  Cd  continuing  to  decrease  slightly  for  higher 
Reynolds  numbers.  Still,  considering  how  much  the  experimental  results  vary 
over  this  range,  our  data  is  not  a  terrible  fit.  Our  model  is  a  reasonable 
approximation  of  dormant  Spartina  alterniflora  above  a  Reynolds  number  of 
450  or  500.  The  flow  depth  appeared  to  have  more  of  an  effect  here  than  in 
the  other  simulations. 

Our  model  did  not  do  a  particularly  good  job  of  capturing  the  drag 
effects  of  live,  green  Spartina  alterniflora.  While  we  capture  the  slight  decrease 
in  Cd  with  increasing  Reynolds  number,  we  greatly  underestimate  its  value  over 
the  entire  range.  The  curve  fit  to  the  SERRI  data  was  calculated  in  the  report 
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Figure  4.21:  Bulk  drag  coefficient  values  for  dormant  Spartina  alterniflora 
vegetation  setup  compared  to  Wu  [115]  results. 

as 

Cd  =  {m67/Ref-\ 

In  the  report  explaining  the  experimental  data,  the  drag  coefficients  are  much 
higher  for  green  Spartina  alterniflora  than  for  any  other  type  of  model  vege¬ 
tation,  and  we  were  unable  to  capture  this. 

It  is  not  surprising  that  our  model  did  not  perfectly  capture  the  drag 
effects  of  tightly  packed  real  vegetation.  Our  model  beds  of  vegetation  are 
imperfect  and  do  not  capture  the  geometry  perfectly.  Also,  there  are  a  great 
deal  of  parameters  that  go  into  the  bending  model,  and  these  all  had  to  be 
estimated  from  empirical  relationships.  Still,  we  were  able  to  model  the  bulk 
drag  coefficient  for  dormant  Spartina  alterniflora  relatively  well  given  these 
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Figure  4.22:  Bulk  drag  coefficient  values  for  green  Spartina  alterniflora  vege¬ 
tation  setup  compared  to  Wu  [115]  results. 

imperfect  modeling  techniques. 

For  the  types  of  model  vegetation  where  we  could  capture  the  geometry 
and  parameters  exactly,  the  rigid  vegetation  and  foam  rubber  flexible  vegeta¬ 
tion,  we  did  an  exceptional  job  at  calculating  the  bulk  drag  coefficients  at 
Reynolds  numbers  above  around  600.  If  a  different  drag  coefficient  for  a  single 
beam  was  used  for  low  Reynolds  numbers,  perhaps,  the  immersed  structure 
method  presented  here  could  capture  the  drag  effects  for  lower  Reynolds  num¬ 
bers  as  well.  The  motion  of  flexible  structures  in  a  wave  tank  is  incredibly 
complex,  and  our  method  has  done  a  good  job  at  capturing  the  motion  and 
its  effect  on  drag. 
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Chapter  5 

Conclusions  and  Discussion 

5.1  Conclusions 

In  the  previous  chapters  we  have  presented  two  methods  for  modeling 
incompressible  flow  through  vegetated  domains.  The  first  method  requires 
one  to  completely  resolve  the  geometry  of  the  vegetation  on  one  hnite  ele¬ 
ment  mesh.  A  high-resolution  locally  conservative  stabilized  hnite  element 
method  is  utilized  to  calculate  the  how  held  through  the  domain  due  to  a 
specihed  hydraulic  gradient.  The  drag  due  to  the  form  of  the  vegetative  ob¬ 
stacles  is  imposed  on  the  how  by  imposing  the  standard  no-slip  condition  on 
the  boundary  with  the  obstacle.  The  corresponding  drag  coefficient  for  a  pe¬ 
riodic  cell  hlled  with  vegetative  obstacles  is  calculated  by  solving  an  LES  cell 
problem  and  using  a  volume  averaging  technique.  The  drag  coefficient  of  the 
vegetated  periodic  cell  is  highly  dependent  on  the  cell  geometry  and  the  how 
velocity.  We  notice  that  the  vegetated  domain  is  a  porous  medium,  similar 
to  those  that  are  commonly  studied  in  the  modeling  of  subsurface  how.  We 
parametrize  the  drag  coefficient  of  the  periodic  vegetated  cell  for  a  large  range 
of  Reynolds  numbers  using  Darcy  and  non-Darcy  upscaling  laws.  Non-Darcy 
laws  are  used  to  upscale  high-velocity  inertial  hows  through  the  porous  do¬ 
main.  We  hnd  that  using  our  upscaling  technique  that  the  common  quadratic 
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Darcy-Forchheimer  law  and  a  general  power  law  are  good  upscaled  models  for 
matching  highly-resolved  LES  simulations.  However,  the  non-Darcy  param¬ 
eters  in  these  upscaling  laws  are  highly  sensitive  to  the  range  of  flows  over 
which  they  are  calculated.  If  one  is  to  use  these  laws  to  parametrize  the  cell 
drag  coefficient,  they  must  make  sure  that  the  range  of  Reynolds  numbers  over 
which  they  expect  their  flow  to  reach  is  included  in  the  high-resolution  LES 
simulations  that  are  performed  to  calculate  the  non-Darcy  parameters. 

The  second  method  that  we  present  for  modeling  flow  through  vege¬ 
tated  regions  is  designed  to  handle  larger-scale  and  much  more  complicated 
domains.  The  fully-resolved  method  above  requires  stem-scale  resolution  of 
the  flow  domain.  This  extremely  high  resolution  greatly  limits  the  scale  of 
domain  of  interest  to  one  that  only  contains  a  few  pieces  of  vegetation.  It 
also  requires  the  vegetative  obstacles  to  be  rigid  which  is  an  unrealistic  as¬ 
sumption  for  many  common  types  of  coastal  vegetation.  Allowing  flexibility 
in  that  model  would  introduce  a  level  of  computational  complexity  that  would 
be  extremely  difficult  to  handle,  even  at  the  small  scale. 

The  immersed  structure  approach  presented  here  works  for  larger-scale, 
more  complicated  vegetated  domains.  While  the  hrst  method  works  on  the 
scale  of  a  few  plants,  this  method  can  be  applied  to  thousands  and  probably 
millions  of  plants.  Rather  than  fully  resolving  the  vegetative  obstacles  on  one 
fluid  hnite  element  mesh  and  imposing  a  no-slip  condition  to  induce  the  drag, 
an  immersed  boundary  approach  is  used.  Separate  fluid  and  structure  meshes 
are  used  in  their  natural  Eulerian  and  Lagrangian  frameworks,  respectively. 
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The  vegetative  obstacles  are  modeled  as  long,  thin  inextensible  flexible  can¬ 
tilever  beams.  Classical  Euler-Bernoulli  beam  is  used  in  the  beam  model.  Two 
numerical  methods  are  developed  for  solving  the  beam  model.  One  is  a  fully 
nonlinear  finite  element  method  that  finds  stable  minima  of  the  beam  energy 
functional  using  Newton’s  method.  The  global  minima  of  this  functional  is  the 
correct  static  state  for  the  beam  system.  The  second  method  is  a  linearized 
finite  element  method  that  uses  Picard  iterations  to  solve  the  Euler-Lagrange 
equations  of  the  beam  energy  functional.  These  equations  are  exactly  the 
equations  for  the  balance  of  moments  and  forces  in  an  elastic  beam.  We  found 
that  combined  with  incremental  loading  that  the  fully  nonlinear  method  is  in¬ 
credibly  robust  for  finding  stable  equilibria  for  beam  systems  under  loads  that 
cause  large  deflections.  Using  between  5  and  10  incremental  loading  steps 
required  the  fewest  total  iterations.  The  linearized  method  is  much  cheaper 
computationally  than  the  nonlinear  method  per  iteration.  For  loads  causing 
relatively  small  deflections,  especially  end  loads,  it  provided  results  just  as  ac¬ 
curate  as  the  nonlinear  method  but  for  a  lower  computational  cost.  However, 
for  loads  causing  moderate  and  large  deflections,  the  linearized  method  has 
a  tendency  to  take  an  extremely  large  number  of  iterations  to  converge  and 
also  to  diverge  or  converge  to  unstable  extrema.  Therefore,  we  found  for  most 
cases  the  fully  nonlinear  beam  solver  to  be  the  best  method  to  use. 

We  couple  the  fluid-structure  system  using  an  immersed  boundary  ap¬ 
proach.  Our  main  objective  was  to  calculate  and  parametrize  the  bulk  drag 
of  many  immersed  structures  on  a  fluid  region.  The  standard  drag  equation 
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is  used  to  calculate  the  drag  force  on  the  fluid  by  the  immersed  structures. 
As  is  the  standard  with  immersed  boundary  methods,  integral  transforms  in¬ 
volving  approximations  to  the  Dirac  delta  are  used  to  map  between  the  fluid 
and  structure  domains.  In  a  method  that  conserves  momentum,  the  momen¬ 
tum  loss  due  to  this  drag  force  enters  the  fluid  model  as  a  local  sink  term. 
The  corresponding  distributed  load  on  the  structure  due  to  the  drag  is  applied 
to  the  beam  model.  A  fully  coupled  steady-state  fluid-structure  interaction 
method  is  formulated.  A  dynamic  method  which  will  be  our  main  tool  is  also 
developed.  This  involves  a  fully  dynamic  LES  flnid  model  and  a  quasi-static 
beam  model. 

The  dynamic  fluid-strncture  interaction  model  was  used  to  model  two 
complex  flow  scenarios.  The  hrst  is  channel  flow  containing  many  rigid  or 
flexible  vegetative  obstacles.  We  calculated  bulk  drag  coefficients  for  these 
flows  that  were  in  the  same  range  as  those  calculated  from  laboratory  experi¬ 
ments.  We  also  captured  trends  that  were  noticed  from  experimental  results. 
These  trends  are  that  the  bulk  drag  coefficient  tends  to  decrease  with  increased 
Reynolds  nnmber  and  with  increased  popnlation  density.  Flexible  vegetative 
obstacles  cause  less  bulk  drag  than  rigid  vegetative  obstacles  of  the  same  di¬ 
mension.  Also,  we  found  that  the  obstacle  thickness  and  length  have  no  major 
effect  on  the  bulk  drag  coefficient  which  was  also  found  by  experiments.  These 
hndings  match  those  from  well-respected  experimental  data. 

The  second  flow  scenario  was  that  of  a  wave  tank.  The  drag  due  to 
vegetation  has  a  large  impact  on  wave  attennation  in  coastal  regions.  Using 
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a  two-phase  flow  method  with  a  conservative  level  set  method  used  to  track 
the  water-air  interface,  we  computationally  modeled  a  wavemaker  that  makes 
linear  plane  waves  and  allows  them  to  propagate  through  a  flume.  Inside 
the  flume,  we  placed  rigid  vegetation,  model  flexible  vegetation  made  of  foam 
rubber  with  known  parameters,  and  approximations  of  real  vegetation.  By  av¬ 
eraging  over  several  wave  periods  and  using  several  wave  amplitudes  and  water 
heights,  we  calculated  bulk  drag  coefficients  for  the  wave  tank  setup  for  large 
ranges  of  Reynolds  numbers.  For  the  rigid  vegetation  and  model  flexible  veg¬ 
etation,  our  calculated  bulk  drag  coefficients  match  up  exceedingly  well  with 
experimental  data.  For  the  approximations  of  real  vegetation,  we  captured  the 
main  trend  of  the  bulk  drag  coefficient  to  decease  with  increasing  Reynolds 
number,  but  were  slightly  off  in  magnitude.  This  discrepancy  is  most  likely 
due  to  the  inaccuracy  of  our  method  at  approximating  a  real  bed  of  vegetation 
and  error  in  the  measurements  that  we  used  to  make  this  approximation. 

5.2  Moving  To  Larger  Scale 

We  have  developed  two  methods  for  simulating  flow  through  vegetated 
regions  at  the  small-scale  and  the  meso-scale.  Both  of  these  methods  provide 
drag  coefficients  that  capture  the  drag  effects  of  a  larger  vegetated  region.  The 
largest  length  scales  where  the  methods  presented  here  could  be  used  to  model 
flow  through  vegetation  directly  is  tens  of  meters.  The  expense  of  resolving 
each  stem  individually,  even  with  the  immersed  structure  method,  becomes 
too  high  at  scales  larger  than  that.  However, the  parametrizations  of  the  drag 
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coefficient  of  flows  through  these  smaller  regions  can  be  applied  to  larger  scale 
simulations. 

The  hrst  method,  which  requires  resolving  all  of  the  geometry  of  the 
vegetation  provides  cell  drag  coefficients  for  each  vegetated  cell.  This  drag 
coefficient  is  a  function  of  the  vegetation  geometry  and  Reynolds  number  and 
can  accurately  be  parameterized  over  a  range  of  Reynolds  number  using  non- 
Darcy  laws.  In  a  large-scale  simulation  the  calculated  drag  coefficient  tensors 
Cd  can  be  used  to  add  a  momentum  sink  term  to  the  large  scale  flow  equations 
modeling  the  momentum  loss  to  drag.  The  sink  term  would  be  of  the  form 

F  =  -^pR|v|Cd(v)v.  (5.1) 

The  sink  term  would  be  applied  to  the  large-scale  momentum  equations  in  the 
region  containing  the  geometry  of  the  cell  for  which  was  calculated.  The 
results  could  be  verihed  by  comparing  to  experimental  results  by  Tsujimoto 
et  ah  [104]. 

The  bulk  drag  coefficient  calculated  using  the  immersed  structure  method 
can  easily  be  applied  to  larger  scale  problems.  It  is  important  to  note  that 
the  bulk  drag  coefficient  is  highly  dependent  on  the  flow  velocity,  the  popu¬ 
lation  density  of  the  vegetation,  and  the  type  of  flow.  We  have  seen  that  the 
bulk  drag  coefficients  for  channel  flow  are  much  smaller  than  those  for  wave 
tank  flows.  We  have  also  presented  two  different  scalings  for  the  bulk  drag 
coefficient,  CdH  and  CdK)  scaled  by  the  flow  depth  and  the  average  vegeta¬ 
tion  height  respectively.  For  3D  large-scale  problems,  the  drag  force  per  unit 
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volume  could  enter  as  a  sink  term  of  the  form 


(5.2) 


everywhere  in  the  horizontal  domain  containing  the  vegetation  that  corre¬ 
sponds  to  CdH-  In  a  perhaps  better  model,  the  drag  due  to  vegetation  would 
only  be  applied  below  the  vegetation  height: 


f 


\Cdh^phN^w\v\  z  <h^ 

0  Z  >  hy. 


(5.3) 


Most  commonly,  for  large-scale  simulations,  2D  depth-averaged  flow 
models  are  used.  One  way  that  the  drag  due  to  vegetation  could  be  modeled 
in  this  scenario  is  with  a  depth  integrated  version  of  (5.2)  where  v  is  replaced  by 
the  depth- averaged  velocity.  As  discussed  by  Westerink  et  ah  [109],  generally 
for  these  large-scale  2D  models,  the  drag  due  to  vegetation  is  incorporated 
into  the  bottom  friction  term 


n  =  Cf[{U^  +  V^y/yH]  (5.4) 

where  r*  is  the  bottom  friction  ,  Cf  is  the  nonlinear  bottom  friction  coefficient, 
H  is  the  water  depth,  and  U  and  V  are  the  depth  averaged  horizontal  velocities. 
Cf  is  commonly  dehned  using  Manning’s  n: 

Cf  =  9^.  (5.5) 

This  idea  was  proposed  by  Ree  and  Palmer  [87].  Usually  the  value  of  Manning’s 
n  is  based  on  extremely  general  land  use  descriptions  and  it  incorporates  both 
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bottom  friction  and  the  form  drag  due  to  vegetation.  Manning’s  equation 
is  the  proper  model  for  bottom  friction,  but  it  is  not  the  proper  model  for 
form  drag.  However,  since  it  often  provides  good  results  for  large-scale  flow 
simulation  the  effect  of  form  drag  due  to  vegetation  is  often  incorporated  in 
this  method  of  ramping  up  the  bottom  friction  parameter  to  incorporate  the 
form  drag.  This  is  discussed  in  depth  by  Wamsley  et  ah  [108]. 

Recently  a  method  has  been  developed  by  Wu  [114]  to  directly  incorpo¬ 
rate  the  bulk  drag  coefficient  due  to  vegetation  that  we  have  discussed  in  this 
dissertation  into  Manning’s  n.  Where  Ub  is  the  bottom  roughness  the  model 

is 

=nl  +  ^  ^  ^^CdKA^ril  min(^,  1.0)/?^^  (5.6) 

where  c^o  is  a  known  parameter,  rj^  is  a  coefficient  of  about  1.0  introduced  by 
Stone  and  Shen  [100],  is  the  frontal  area  of  the  vegetation,  and  Rg  is  the 
hydraulic  gradient.  (5.6)  builds  off  of  a  model  originally  proposed  by  Petryk 
and  Bosmajian  [79]  where  they  modeled  flow  through  mangroves  by  increasing 
Manning’s  n  based  on  stem  population  densities  and  thicknesses.  Other  work 
using  these  types  of  Manning’s  n  models  incorporating  the  bulk  drag  effects 
with  mangroves  has  been  done  by  Wolanski  et  al.  [113],  Mazda  et  al.  [66],  and 
Furukawa  et  al.  [28].  Drag  coefficients  calculated  using  the  methods  presented 
here  could  be  incorporated  into  a  Manning’s  n  model  of  this  type  relatively 
easily  allowing  it  to  be  used  for  large-scale  2D  simulations  like  those  done  using 
the  Advanced  Circulation  (ADCIRC)  model  [63]. 
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5.3  Further  Work 

There  is  a  great  deal  of  future  work  to  be  done  in  the  modeling  of 
flow  through  vegetative  regions.  Incorporating  the  drag  coefficient  calcula¬ 
tions  presented  in  this  dissertation  into  large-scale  2D  and  3D  simulations 
as  discussed  in  the  previous  section  is  an  obvious  direction  for  future  work. 
These  large-scale  methods  can  be  used  to  model  important  flows  in  riverine 
and  coastal  regions.  These  areas  are  of  utmost  importance  for  a  variety  of 
environmental,  industrial,  economic,  and  defense  applications.  In  these  areas, 
small-scale  effects,  including  the  effect  of  drag  due  to  vegetation  play  a  large 
role  in  determining  large-scale  flow  characteristics.  If  the  methods  for  quan¬ 
tifying  vegetative  drag  presented  here,  which  we  have  shown  to  be  accurate, 
can  improve  larger-scale  models  it  would  be  an  important  contribution  to  the 
held  of  hydrology. 

In  this  dissertation,  we  have  mostly  discussed  vegetation  interacting 
with  incompressible  hows  through  drag  ehects.  Drag  surely  is  the  ehect  that 
has  the  greatest  ehect  on  mean  how;  however,  vegetation  also  greatly  impacts 
other  how  processes  including  wave  attenuation,  turbulence,  dihusion,  and 
transport.  We  have  already  discussed  the  ehect  of  vegetation  on  the  atten¬ 
uation  of  waves,  but  only  with  simple  linear  waves.  There  is  a  great  deal  of 
possible  work  to  be  done  studying  the  ehect  of  vegetation  on  breaking  waves 
and  other  more  complicated  wave  phenomena.  Most  work  in  this  area  has  been 
theoretical  or  experimental,  but  we  have  shown  that  computational  methods 
can  be  used  to  model  complex  physical  processes  such  as  these. 
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The  existence  of  vegetation  in  a  high-velocity  flow  greatly  increases  the 
production  of  turbulence.  We  have  utilized  a  large  eddy  simulation  flow  model 
in  our  calculations  in  order  to  capture  some  of  the  effects  due  to  this  wake 
production.  Nepf  et  al.  [74]  and  Burke  and  Stolzenbach  [11],  show  that  even 
for  sparsely  populated  vegetation,  turbulence  production  due  to  vegetation  is 
greater  than  that  due  to  bed  shear.  Assuming  that  the  energy  leaving  the 
mean  flow  due  to  drag  becomes  turbulent  kinetic  energy,  then  the  production 
of  turbulent  kinetic  energy  is 

P  =  (5.7) 

According  to  Tennekes  and  Lumley  [101],  the  corresponding  dissipation  of 
turbulent  kinetic  energy  scales  like 

e  ~  (5.8) 

where  k  is  the  turbulent  kinetic  energy.  Using  these  observations,  our  calcu¬ 
lations  of  Cd  could  be  highly  useful  in  modeling  turbulence  due  to  vegetation. 
It  could  also  be  seen  if  using  our  resolved  pore-scale  and  immersed  structure 
methods  are  successful  in  capturing  turbulent  stress  prohles  in  the  flow.  These 
types  of  studies  have  been  done  in  flume  experiments  by  Nepf  and  Ghisalberti 
[72],  Wilson  et  al.  [110],  Dunn  et  al.  [23],  and  Jarvela  [41]  and  computational 
experiments  by  Li  and  Xie  [54]  and  Erduran  and  Kutija  [25]. 

Within  a  flow  through  a  vegetated  domain,  diffusion  can  be  separated 
into  turbulent  diffusion  and  mechanical  diffusion.  Turbulent  diffusion  is  dif¬ 
fusion  caused  by  turbulence.  As  discussed  previously,  vegetation  in  a  flow 
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greatly  increases  the  production  of  turbulence,  so  it  also  greatly  affects  turbu¬ 
lent  diffusion.  Mechanical  diffusion  is  commonly  studied  with  respect  to  flow 
through  porous  media.  According  to  Nepf  [73],  it  concerns  the  dispersal  of 
fluid  particles  because  of  the  complex  flow  paths  due  to  the  porous  nature  of 
the  domain.  As  we  have  shown,  a  vegetated  region  is  a  type  of  porous  struc¬ 
ture,  so  mechanical  diffusion  is  a  concern  in  these  regions.  Our  models  could 
likely  be  applied  to  study  these  diffusion  processes. 

When  studying  flow  processes  it  is  ubiquitous  to  also  study  transport 
phenomena  that  occur  within  that  flow.  Transport  through  a  vegetated  re¬ 
gion  is  an  extremely  complex  process  the  study  of  is  one  for  which  there  is 
room  for  a  great  deal  of  new  research.  Models  for  transport  through  vege¬ 
tated  regions  could  be  used  to  model  contaminants,  nutrients,  fish  eggs,  and 
several  other  things  of  environmental,  ecological,  and  industrial  importance. 
Besides  vegetated  coastal  regions,  the  numerical  modeling  of  flow  and  trans¬ 
port  through  obstructed  flows  may  be  applicable  to  several  other  engineering 
concerns  including  flow  and  transport  through  urban  environments,  forests, 
fields  of  crops,  biohlm  reactors,  and  porous  media. 
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